Index: DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/SoilProfile2D.cs =================================================================== diff -u -r5251 -r5252 --- DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/SoilProfile2D.cs (.../SoilProfile2D.cs) (revision 5251) +++ DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/SoilProfile2D.cs (.../SoilProfile2D.cs) (revision 5252) @@ -71,6 +71,8 @@ /// /// Gets the soil profile 1D at the given X. + /// If the given 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 /// /// The x. /// Soil Profile 1D @@ -94,14 +96,12 @@ { x = Geometry.Left + diff; } - - // 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 { - isXShiftedToRight = HasLayerAVerticalPartAtGivenX(xRelevant) && !AreSoilProfilesIdentical(xRelevant - diff / 2, xRelevant + diff / 2); + isXShiftedToRight = AreLayersWithVerticalPartPresentAtGivenX(xRelevant) && !AreSoilProfilesIdentical(xRelevant - diff / 2, xRelevant + diff / 2); if (isXShiftedToRight) { xRelevant += diff; @@ -156,7 +156,7 @@ public static SoilLayer2D GetOriginalLayerFromOldSurfaces(GeometrySurface geometrySurface, IEnumerable oldSurfaces) { - return GetOriginalLayerFromOldSurfaces(geometrySurface, oldSurfaces, 0.0); + return GetOriginalLayerFromOldSurfaces(geometrySurface, oldSurfaces, 0.0); } /// @@ -169,29 +169,29 @@ private static SoilLayer2D GetOriginalLayerFromOldSurfaces(GeometrySurface geometrySurface, IEnumerable oldSurfaces, double shift) { - Point2D point = new Point2D(0.0, 0.0); - bool isPointInOuterLoopAndOldSurface = false; - point = IsPointInOuterLoopAndOldSurface(geometrySurface, oldSurfaces, shift, point, ref isPointInOuterLoopAndOldSurface); + Point2D point = new Point2D(0.0, 0.0); + bool isPointInOuterLoopAndOldSurface = false; + point = IsPointInOuterLoopAndOldSurface(geometrySurface, oldSurfaces, shift, point, ref isPointInOuterLoopAndOldSurface); - if (!isPointInOuterLoopAndOldSurface) - { - isPointInOuterLoopAndOldSurface = IsPointInOldSurfaceJustBelowTopOfNewGeometryWithinItsLimits(geometrySurface, - oldSurfaces, shift, point, isPointInOuterLoopAndOldSurface); - } - if (isPointInOuterLoopAndOldSurface) - { - point.X -= shift; - if (IsPointInPreviousOuterLoopOfOldSurface(oldSurfaces, point, out SoilLayer2D soilLayer2D)) + if (!isPointInOuterLoopAndOldSurface) { - return soilLayer2D; + isPointInOuterLoopAndOldSurface = IsPointInOldSurfaceJustBelowTopOfNewGeometryWithinItsLimits(geometrySurface, + oldSurfaces, shift, point, isPointInOuterLoopAndOldSurface); } - - if (IsPointInOuterLoopOfOldSurface(oldSurfaces, point, out SoilLayer2D originalLayerFromOldSurfaces1)) + if (isPointInOuterLoopAndOldSurface) { - return originalLayerFromOldSurfaces1; + point.X -= shift; + if (IsPointInPreviousOuterLoopOfOldSurface(oldSurfaces, point, out SoilLayer2D soilLayer2D)) + { + return soilLayer2D; + } + + if (IsPointInOuterLoopOfOldSurface(oldSurfaces, point, out SoilLayer2D originalLayerFromOldSurfaces1)) + { + return originalLayerFromOldSurfaces1; + } } - } - return null; + return null; } private static bool IsPointInOuterLoopOfOldSurface(IEnumerable oldSurfaces, Point2D point, @@ -202,7 +202,7 @@ { GeometryLoop outerLoop = oldSurface.GeometrySurface.OuterLoop; if (outerLoop != null && outerLoop.CurveList.Count > 2 && Routines2D.CheckIfPointIsInPolygon(outerLoop, point.X, - point.Z) == PointInPolygon.InsidePolygon) + point.Z) == PointInPolygon.InsidePolygon) { bool isPointInOuterLoopOfOldSurface = true; foreach (GeometryLoop innerLoop in oldSurface.GeometrySurface.InnerLoops) @@ -249,7 +249,7 @@ } private static bool IsPointInOldSurfaceJustBelowTopOfNewGeometryWithinItsLimits(GeometrySurface geometrySurface, - IEnumerable oldSurfaces, double shift, Point2D point, bool isPointInOuterLoopAndOldSurface) + IEnumerable oldSurfaces, double shift, Point2D point, bool isPointInOuterLoopAndOldSurface) { GeometryPointString topGeometrySurface = geometrySurface.DetermineTopGeometrySurface(); topGeometrySurface.SortPointsByXAscending(); @@ -365,7 +365,7 @@ foreach (GeometryLoop innerLoop in innerLoops) { if (Routines2D.CheckIfPointIsInPolygon(innerLoop, point.X, point.Z) == PointInPolygon.InsidePolygon) - return false; + return false; } } return isPointInSurface; @@ -400,21 +400,11 @@ return soilProfile; } - private bool HasLayerAVerticalPartAtGivenX(double xCoordinate) + private bool AreLayersWithVerticalPartPresentAtGivenX(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; + return Surfaces.Any(surface => surface.GeometrySurface.OuterLoop.CurveList.Any + (curve => curve.HeadPoint.X.IsNearEqual(xCoordinate, toleranceAlmostEqual) && + curve.EndPoint.X.IsNearEqual(xCoordinate, toleranceAlmostEqual))); } private bool AreSoilProfilesIdentical(double x1, double x2) Index: DamEngine/trunk/src/Deltares.DamEngine.Calculators.Tests/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverterTests.cs =================================================================== diff -u -r5120 -r5252 --- DamEngine/trunk/src/Deltares.DamEngine.Calculators.Tests/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverterTests.cs (.../PlLinesToWaternetConverterTests.cs) (revision 5120) +++ DamEngine/trunk/src/Deltares.DamEngine.Calculators.Tests/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverterTests.cs (.../PlLinesToWaternetConverterTests.cs) (revision 5252) @@ -42,7 +42,7 @@ private const string mapTestFiles = @"KernelWrappers\MacroStabilityCommon\TestFiles\"; private const double precision3Decimals = 0.00051; - private const double precision6Decimals = 0.00000051; + private const double precision5Decimals = 0.0000051; private const double penetrationLength = 0.5; private static readonly Random random = new(21); @@ -1859,8 +1859,8 @@ GeometryPoint actualPoint = actualPoints.ElementAt(i); Assert.Multiple(() => { - Assert.That(actualPoint.X, Is.EqualTo(expectedPoint.X).Within(precision6Decimals)); - Assert.That(actualPoint.Z, Is.EqualTo(expectedPoint.Z).Within(precision6Decimals)); + Assert.That(actualPoint.X, Is.EqualTo(expectedPoint.X).Within(precision5Decimals)); + Assert.That(actualPoint.Z, Is.EqualTo(expectedPoint.Z).Within(precision5Decimals)); }); } } Index: DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs =================================================================== diff -u -r5085 -r5252 --- DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs (.../PlLinesToWaternetConverter.cs) (revision 5085) +++ DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs (.../PlLinesToWaternetConverter.cs) (revision 5252) @@ -43,7 +43,9 @@ LowestLayer } - private const double epsilon = 1e-9; + private const double deviationX = 1e-06; + private const double toleranceAlmostEqual = 1e-09; + private const string waternetLine1Name = "Waternet line phreatic line"; private const string waternetLine2Name = "Penetration zone lower aquifer"; private const string waternetLine3Name = "Waternet line lower aquifer"; @@ -333,7 +335,7 @@ } line.SyncCalcPoints(); - LineHelper.RemoveDuplicatedPoints(line.CalcPoints, epsilon); + LineHelper.RemoveDuplicatedPoints(line.CalcPoints, toleranceAlmostEqual); line.RemoveUnnecessaryPoints(); line.SyncPoints(); @@ -361,8 +363,8 @@ { if (HasAquiferAVerticalPartAtGivenX(layerType, xCoordinate, soilProfile, indexInBetweenAquifer)) { - xCoordinatesAll.Add(xCoordinate - epsilon); - xCoordinatesAll.Add(xCoordinate + epsilon); + xCoordinatesAll.Add(xCoordinate - deviationX); + xCoordinatesAll.Add(xCoordinate + deviationX); } else { @@ -403,8 +405,8 @@ private static bool IsLayerBoundaryContinuous(List boundaryPoints, double leftGeometryBoundary, double rightGeometryBoundary) { - return boundaryPoints.Any() && (Math.Abs(boundaryPoints.First().X - leftGeometryBoundary) < epsilon) - && Math.Abs(boundaryPoints.Last().X - rightGeometryBoundary) < epsilon; + return boundaryPoints.Any() && (Math.Abs(boundaryPoints.First().X - leftGeometryBoundary) < toleranceAlmostEqual) + && Math.Abs(boundaryPoints.Last().X - rightGeometryBoundary) < toleranceAlmostEqual; } private static SoilLayer1D GetSoilLayer1D(LayerType layerType, SoilProfile1D soilProfile1D, int indexInBetweenAquifer) @@ -565,7 +567,7 @@ } // For lowest layers that intersect the phreatic line the bottom of these layers is kept - if (intersectionPoints.Count > 0 && outerLoop.Points.Any(p => Math.Abs(p.Z - soilProfile2D.GetSoilProfile1D(p.X).BottomLevel) < epsilon)) + if (intersectionPoints.Count > 0 && outerLoop.Points.Any(p => Math.Abs(p.Z - soilProfile2D.GetSoilProfile1D(p.X).BottomLevel) < toleranceAlmostEqual)) { var geometrySurface = new GeometrySurface(outerLoop); GeometryPointString bottom = geometrySurface.DetermineBottomGeometrySurface(); @@ -582,11 +584,11 @@ List splitLines = SplitLinesAtXCoordinates(xCoordinates, lines); for (var i = 0; i < xCoordinates.Length - 1; i++) { - List relevantLines = splitLines.FindAll(line => Math.Abs(line.BeginPoint.X - xCoordinates[i]) < epsilon && Math.Abs(line.BeginPoint.X - line.EndPoint.X) > 0); + List relevantLines = splitLines.FindAll(line => Math.Abs(line.BeginPoint.X - xCoordinates[i]) < toleranceAlmostEqual && Math.Abs(line.BeginPoint.X - line.EndPoint.X) > 0); double maxZBegin = relevantLines.Select(line => line.BeginPoint.Z).Max(); - List relevantLines2 = relevantLines.FindAll(line => Math.Abs(line.BeginPoint.Z - maxZBegin) < epsilon); + List relevantLines2 = relevantLines.FindAll(line => Math.Abs(line.BeginPoint.Z - maxZBegin) < toleranceAlmostEqual); double maxZEnd = relevantLines2.Select(line => line.EndPoint.Z).Max(); - Line relevantLine = relevantLines2.Find(line => Math.Abs(line.EndPoint.Z - maxZEnd) < epsilon); + Line relevantLine = relevantLines2.Find(line => Math.Abs(line.EndPoint.Z - maxZEnd) < toleranceAlmostEqual); waternetLine.Points.Add(new GeometryPoint(relevantLine.BeginPoint.X, relevantLine.BeginPoint.Z)); waternetLine.Points.Add(new GeometryPoint(relevantLine.EndPoint.X, relevantLine.EndPoint.Z)); } @@ -599,7 +601,7 @@ } waternetLine.SyncCalcPoints(); - LineHelper.RemoveDuplicatedPoints(waternetLine.CalcPoints, epsilon); + LineHelper.RemoveDuplicatedPoints(waternetLine.CalcPoints, toleranceAlmostEqual); waternetLine.RemoveUnnecessaryPoints(); waternetLine.SyncPoints(); @@ -608,8 +610,8 @@ private static bool HasAquiferAVerticalPartAtGivenX(LayerType layerType, double xCoordinate, SoilProfile2D soilProfile, int indexInBetweenAquifer) { - SoilProfile1D crossSectionLeft = soilProfile.GetSoilProfile1D(xCoordinate - epsilon); - SoilProfile1D crossSectionRight = soilProfile.GetSoilProfile1D(xCoordinate + epsilon); + SoilProfile1D crossSectionLeft = soilProfile.GetSoilProfile1D(xCoordinate - deviationX); + SoilProfile1D crossSectionRight = soilProfile.GetSoilProfile1D(xCoordinate + deviationX); if (crossSectionLeft != null && crossSectionRight != null) { SoilLayer1D currentAquiferLeft = GetSoilLayer1D(layerType, crossSectionLeft, indexInBetweenAquifer); @@ -714,7 +716,7 @@ foreach (GeometryPoint waternetLine1Point in from waternetLineOther in waternet.WaternetLineList.FindAll(w => w.Name != waternetLine1Name) from waternetLine1Point in waternetLine1.Points - where waternetLineOther.Points.Any(point => Math.Abs(point.Z - waternetLine1Point.Z) < epsilon) + where waternetLineOther.Points.Any(point => Math.Abs(point.Z - waternetLine1Point.Z) < toleranceAlmostEqual) select waternetLine1Point) { waternetLine1Point.Z += minimumDistance; @@ -728,7 +730,7 @@ var splitLines = new List(); foreach (Line line in lines) { - List xCoordinatesToBeAdded = xCoordinates.Where(xCoordinate => line.BeginPoint.X.IsLessThan(xCoordinate, epsilon) && xCoordinate.IsLessThan(line.EndPoint.X, epsilon)).ToList(); + List xCoordinatesToBeAdded = xCoordinates.Where(xCoordinate => line.BeginPoint.X.IsLessThan(xCoordinate, toleranceAlmostEqual) && xCoordinate.IsLessThan(line.EndPoint.X, toleranceAlmostEqual)).ToList(); if (xCoordinatesToBeAdded.Count > 0) {