Index: DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs =================================================================== diff -u -r6245 -r6404 --- DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs (.../PlLinesToWaternetConverter.cs) (revision 6245) +++ DamEngine/trunk/src/Deltares.DamEngine.Calculators/KernelWrappers/MacroStabilityCommon/PlLinesToWaternetConverter.cs (.../PlLinesToWaternetConverter.cs) (revision 6404) @@ -1,4 +1,4 @@ -// Copyright (C) Stichting Deltares 2024. All rights reserved. +// Copyright (C) Stichting Deltares 2025. All rights reserved. // // This file is part of the Dam Engine. // @@ -31,14 +31,13 @@ using Deltares.DamEngine.Data.Geotechnics; using Deltares.DamEngine.Data.Standard; using CharacteristicPointType = Deltares.DamEngine.Data.Geotechnics.CharacteristicPointType; -using SoilProfile = Deltares.DamEngine.Data.Geotechnics.SoilProfile; namespace Deltares.DamEngine.Calculators.KernelWrappers.MacroStabilityCommon; public static class PlLinesToWaternetConverter { 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"; @@ -66,7 +65,7 @@ // The validation of the soil profile was already performed in the Pl-lines creator. // Here we only have to check that the PL-lines were indeed created. ThrowWhenPlLinesIsNull(plLines); - + double xLeft = surfaceLine.CharacteristicPoints.GetPoint2D(CharacteristicPointType.SurfaceLevelOutside).X; double xRight = surfaceLine.CharacteristicPoints.GetPoint2D(CharacteristicPointType.SurfaceLevelInside).X; @@ -86,7 +85,7 @@ double x = surfaceLine.CharacteristicPoints.GetPoint2D(CharacteristicPointType.DikeToeAtPolder).X; var dummySoil = new Soil("dummy", 0.0, 0.0); SoilProfile1D soilProfile1DDikeToeAtPolder = SoilProfileHelper.DetermineForSurfaceLineCorrected1DProfileAtX(soilProfile1D, surfaceLine, x, dummySoil); - + plLine = plLines.Lines[PlLineType.Pl2]; headLine = CreateLine(plLine, headLine2Name); if (headLine != null && plLine.Points.Count != 0 && @@ -136,9 +135,9 @@ // The validation of the soil profile was already performed in the Pl-lines creator. // Here we only have to check that the PL-lines were indeed created. ThrowWhenPlLinesIsNull(plLines); - + SoilProfile2DHelper.DetermineAquiferLayerBoundaryPoints(SoilProfile2DHelper.LayerType.BottomAquiferCluster, soilProfile, out Point2D[] bottomAquiferCoordinates, out _); - + var waternet = new Waternet(); PlLine plLine = plLines.Lines[PlLineType.Pl1]; if (plLine != null && plLine.Points.Count != 0) @@ -188,15 +187,132 @@ waternetLine.Points.Add(new Point2D(xRight, level)); return waternetLine; } - + + /// + /// Create the waternet line associated to the phreatic line (PL1) for the waternet type "Dam Standard". + /// This line lies on the bottom of the soil layers “in which the phreatic plane lies”. + /// Where PL1 lies above the ground level, the line lies on the ground level. + /// As a consequence, the line can contain vertical "jumps" when layers have vertical parts. + /// + /// + /// + /// + internal static WaternetLine CreateWaternetLineForPhreaticLineForDamStandard(SoilProfile2D soilProfile2D, PlLine plLine) + { + var waternetLine = new WaternetLine(); + var curves = new List(); + foreach (SoilLayer2D surface in soilProfile2D.Surfaces) + { + GeometryLoop outerLoop = surface.GeometrySurface.OuterLoop; + // Add the first point to get a closed GeometryPointString + outerLoop.Points.Add(new Point2D(outerLoop.Points[0].X, outerLoop.Points[0].Z)); + + var plLineInPoint2D = new GeometryPointString(); + foreach (PlLinePoint point in plLine.Points) + { + plLineInPoint2D.Points.Add(new Point2D(point.X, point.Z)); + } + + List intersectionPoints = outerLoop.IntersectionXzPointsWithGeometryString(plLineInPoint2D); + // Only the layers that intersect the phreatic line are kept. + if (intersectionPoints.Count > 0) + { + var geometrySurface = new GeometrySurface(outerLoop); + GeometryPointString bottom = geometrySurface.DetermineBottomGeometrySurface(); + bottom.SortPointsByXAscending(); + curves.Add(bottom); + } + + // Add also the surface line in case of free water + curves.Add(soilProfile2D.Geometry.SurfaceLine); + } + + // The waternet line is the polyline formed with the lowest lines connected to each other. + List lines = DivideCurvesIntoLines(curves); + IEnumerable allPoints = curves.SelectMany(curve => curve.Points); + double[] xCoordinates = allPoints.Select(point => point.X).OrderBy(x => x).Distinct().ToArray(); + List splitLines = SplitLinesAtXCoordinates(xCoordinates, lines); + for (var i = 0; i < xCoordinates.Length - 1; i++) + { + List relevantLinesAtX = splitLines.FindAll(line => Math.Abs(line.BeginPoint.X - xCoordinates[i]) < toleranceAlmostEqual && Math.Abs(line.BeginPoint.X - line.EndPoint.X) > 0); + if (relevantLinesAtX.Count == 0) + { + continue; + } + + double minZBegin = relevantLinesAtX.Select(line => line.BeginPoint.Z).Min(); + List relevantLinesAtXWithZBeginMax = relevantLinesAtX.FindAll(line => Math.Abs(line.BeginPoint.Z - minZBegin) < toleranceAlmostEqual); + double minZEnd = relevantLinesAtXWithZBeginMax.Select(line => line.EndPoint.Z).Min(); + Line relevantLine = relevantLinesAtXWithZBeginMax.Find(line => Math.Abs(line.EndPoint.Z - minZEnd) < toleranceAlmostEqual); + waternetLine.Points.Add(new Point2D(relevantLine.BeginPoint.X, relevantLine.BeginPoint.Z)); + waternetLine.Points.Add(new Point2D(relevantLine.EndPoint.X, relevantLine.EndPoint.Z)); + } + + // Make sure the waternet line for PL 1 is not in the bottom aquifer. + foreach (Point2D waternetLinePoint in waternetLine.Points) + { + double bottomAquiferLevel = soilProfile2D.GetSoilProfile1D(waternetLinePoint.X).BottomAquiferLayer.TopLevel; + waternetLinePoint.Z = Math.Max(waternetLinePoint.Z, bottomAquiferLevel); + } + + LineHelper.RemoveDuplicatedPoints(waternetLine.Points, toleranceAlmostEqual); + waternetLine.RemoveUnnecessaryPoints(); + AvoidUpwardsVerticalLine(waternetLine.Points); + + return waternetLine; + } + + internal static void SplitLineAtXCoordinate(List xCoordinatesToBeAdded, Line line, List splitLines) + { + for (var i = 0; i < xCoordinatesToBeAdded.Count; i++) + { + var point1 = new Point2D(xCoordinatesToBeAdded[i], Math.Max(line.BeginPoint.Z, line.EndPoint.Z) + 1); + var point2 = new Point2D(xCoordinatesToBeAdded[i], Math.Min(line.BeginPoint.Z, line.EndPoint.Z) - 1); + double zCoordinate = line.GetIntersectPointXz(new Line(point1, point2)).Z; + + var middlePoint = new Point2D(xCoordinatesToBeAdded[i], zCoordinate); + Point2D beginPoint = i == 0 ? line.BeginPoint : splitLines.Last().EndPoint; + + splitLines.Add(new Line(beginPoint, middlePoint)); + if (i == xCoordinatesToBeAdded.Count - 1) + { + splitLines.Add(new Line(middlePoint, line.EndPoint)); + } + } + } + + /// + /// This correction is needed for the creation of the STIX file because a reference line containing a vertical part + /// is always drawn in downwards direction by D-Stability. + /// + internal static void AvoidUpwardsVerticalLine(IList points) + { + for (var i = 0; i < points.Count - 1; i++) + { + if (!points[i].X.IsNearEqual(points[i + 1].X, toleranceAlmostEqual) || points[i].Z.IsGreaterThanOrEqualTo(points[i + 1].Z, toleranceAlmostEqual)) + { + continue; + } + + if (i == points.Count - 2 || (points[i + 1].X + 0.01).IsLessThan(points[i + 2].X, toleranceAlmostEqual)) + { + points[i + 1].X += 0.01; + } + else + { + points[i + 1].X = points[i + 2].X - (points[i + 2].X - points[i + 1].X) / 2; + } + } + } + private static void CreateWaternetLinesForInBetweenAquifers(Waternet waternet, SoilProfile1D soilProfile1D, double xLeft, double xRight, PlLine plLine4, IntrusionVerticalWaterPressureType? pressureType) { if (plLine4 == null || plLine4.Points.Count == 0 || soilProfile1D.GetInBetweenAquiferClusters == null || pressureType == IntrusionVerticalWaterPressureType.FullHydroStatic) { return; } - + var headLine = CreateLine(plLine4, headLine4Name); waternet.HeadLineList.Add(headLine); var number = 0; @@ -208,7 +324,7 @@ waternetLineTop.Name = waternetLine4TopName + " (" + number + ")"; waternetLineTop.HeadLine = headLine; waternet.WaternetLineList.Add(waternetLineTop); - + double levelBottom = inBetweenAquiferCluster.Item2.BottomLevel; WaternetLine waternetLineBottom = CreateWaternetLine(levelBottom, xLeft, xRight); waternetLineBottom.Name = waternetLine4BottomName + " (" + number + ")"; @@ -225,7 +341,7 @@ return; } - SoilProfile2DHelper.DetermineInBetweenAquifersLayerBoundaryPoints(soilProfile, out List inBetweenAquiferUpperCoordinates, + SoilProfile2DHelper.DetermineInBetweenAquifersLayerBoundaryPoints(soilProfile, out List inBetweenAquiferUpperCoordinates, out List inBetweenAquiferLowerCoordinates); if (!inBetweenAquiferUpperCoordinates.Any() || !inBetweenAquiferLowerCoordinates.Any()) @@ -282,10 +398,10 @@ LineHelper.RemoveDuplicatedPoints(line.Points, toleranceAlmostEqual); line.RemoveUnnecessaryPoints(); AvoidUpwardsVerticalLine(line.Points); - + return line; } - + /// /// The position of the waternet line associated to the phreatic line depends on the water pressure type: /// - For DAM Standard type, this line lies on the bottom of the soil layers “in which the phreatic plane lies”. @@ -391,79 +507,6 @@ } /// - /// Create the waternet line associated to the phreatic line (PL1) for the waternet type "Dam Standard". - /// This line lies on the bottom of the soil layers “in which the phreatic plane lies”. - /// Where PL1 lies above the ground level, the line lies on the ground level. - /// As a consequence, the line can contain vertical "jumps" when layers have vertical parts. - /// - /// - /// - /// - internal static WaternetLine CreateWaternetLineForPhreaticLineForDamStandard(SoilProfile2D soilProfile2D, PlLine plLine) - { - var waternetLine = new WaternetLine(); - var curves = new List(); - foreach (SoilLayer2D surface in soilProfile2D.Surfaces) - { - GeometryLoop outerLoop = surface.GeometrySurface.OuterLoop; - // Add the first point to get a closed GeometryPointString - outerLoop.Points.Add(new Point2D(outerLoop.Points[0].X, outerLoop.Points[0].Z)); - - var plLineInPoint2D = new GeometryPointString(); - foreach (PlLinePoint point in plLine.Points) - { - plLineInPoint2D.Points.Add(new Point2D(point.X, point.Z)); - } - - List intersectionPoints = outerLoop.IntersectionXzPointsWithGeometryString(plLineInPoint2D); - // Only the layers that intersect the phreatic line are kept. - if (intersectionPoints.Count > 0) - { - var geometrySurface = new GeometrySurface(outerLoop); - GeometryPointString bottom = geometrySurface.DetermineBottomGeometrySurface(); - bottom.SortPointsByXAscending(); - curves.Add(bottom); - } - - // Add also the surface line in case of free water - curves.Add(soilProfile2D.Geometry.SurfaceLine); - } - - // The waternet line is the polyline formed with the lowest lines connected to each other. - List lines = DivideCurvesIntoLines(curves); - IEnumerable allPoints = curves.SelectMany(curve => curve.Points); - double[] xCoordinates = allPoints.Select(point => point.X).OrderBy(x => x).Distinct().ToArray(); - List splitLines = SplitLinesAtXCoordinates(xCoordinates, lines); - for (var i = 0; i < xCoordinates.Length - 1; i++) - { - List relevantLinesAtX = splitLines.FindAll(line => Math.Abs(line.BeginPoint.X - xCoordinates[i]) < toleranceAlmostEqual && Math.Abs(line.BeginPoint.X - line.EndPoint.X) > 0); - if (relevantLinesAtX.Count == 0) - { - continue; - } - double minZBegin = relevantLinesAtX.Select(line => line.BeginPoint.Z).Min(); - List relevantLinesAtXWithZBeginMax = relevantLinesAtX.FindAll(line => Math.Abs(line.BeginPoint.Z - minZBegin) < toleranceAlmostEqual); - double minZEnd = relevantLinesAtXWithZBeginMax.Select(line => line.EndPoint.Z).Min(); - Line relevantLine = relevantLinesAtXWithZBeginMax.Find(line => Math.Abs(line.EndPoint.Z - minZEnd) < toleranceAlmostEqual); - waternetLine.Points.Add(new Point2D(relevantLine.BeginPoint.X, relevantLine.BeginPoint.Z)); - waternetLine.Points.Add(new Point2D(relevantLine.EndPoint.X, relevantLine.EndPoint.Z)); - } - - // Make sure the waternet line for PL 1 is not in the bottom aquifer. - foreach (Point2D waternetLinePoint in waternetLine.Points) - { - double bottomAquiferLevel = soilProfile2D.GetSoilProfile1D(waternetLinePoint.X).BottomAquiferLayer.TopLevel; - waternetLinePoint.Z = Math.Max(waternetLinePoint.Z, bottomAquiferLevel); - } - - LineHelper.RemoveDuplicatedPoints(waternetLine.Points, toleranceAlmostEqual); - waternetLine.RemoveUnnecessaryPoints(); - AvoidUpwardsVerticalLine(waternetLine.Points); - - return waternetLine; - } - - /// /// Throws when the pl lines object is not assigned. /// /// The pl lines. @@ -505,7 +548,7 @@ foreach (Line line in lines) { List xCoordinatesToBeAdded = xCoordinates.Where(xCoordinate => line.BeginPoint.X.IsLessThan(xCoordinate, toleranceAlmostEqual) && xCoordinate.IsLessThan(line.EndPoint.X, toleranceAlmostEqual)).ToList(); - + if (xCoordinatesToBeAdded.Count > 0) { SplitLineAtXCoordinate(xCoordinatesToBeAdded, line, splitLines); @@ -519,25 +562,6 @@ return splitLines; } - internal static void SplitLineAtXCoordinate(List xCoordinatesToBeAdded, Line line, List splitLines) - { - for (var i = 0; i < xCoordinatesToBeAdded.Count; i++) - { - var point1 = new Point2D(xCoordinatesToBeAdded[i], Math.Max(line.BeginPoint.Z, line.EndPoint.Z) + 1); - var point2 = new Point2D(xCoordinatesToBeAdded[i], Math.Min(line.BeginPoint.Z, line.EndPoint.Z) - 1); - double zCoordinate = line.GetIntersectPointXz(new Line(point1, point2)).Z; - - var middlePoint = new Point2D(xCoordinatesToBeAdded[i], zCoordinate); - Point2D beginPoint = i == 0 ? line.BeginPoint : splitLines.Last().EndPoint; - - splitLines.Add(new Line(beginPoint, middlePoint)); - if (i == xCoordinatesToBeAdded.Count - 1) - { - splitLines.Add(new Line(middlePoint, line.EndPoint)); - } - } - } - private static List DivideCurvesIntoLines(List curves) { var lines = new List(); @@ -557,30 +581,6 @@ return lines; } - /// - /// This correction is needed for the creation of the STIX file because a reference line containing a vertical part - /// is always drawn in downwards direction by D-Stability. - /// - internal static void AvoidUpwardsVerticalLine(IList points) - { - for (var i = 0; i < points.Count - 1; i++) - { - if (!points[i].X.IsNearEqual(points[i + 1].X, toleranceAlmostEqual) || points[i].Z.IsGreaterThanOrEqualTo(points[i + 1].Z, toleranceAlmostEqual)) - { - continue; - } - - if (i == points.Count - 2 || (points[i + 1].X + 0.01).IsLessThan(points[i + 2].X, toleranceAlmostEqual)) - { - points[i + 1].X += 0.01; - } - else - { - points[i + 1].X = points[i + 2].X - (points[i + 2].X - points[i + 1].X) / 2; - } - } - } - private static bool IsLineEndingInPolderSide(GeometryPointString line, SurfaceLine2 surfaceLine) { Point2D dikeToeAtPolder = surfaceLine.CharacteristicPoints.GetPoint2D(CharacteristicPointType.DikeToeAtPolder);