// Copyright (C) Stichting Deltares 2024. All rights reserved. // // This file is part of the Dam Engine. // // The Dam Engine is free software: you can redistribute it and/or modify // it under the terms of the GNU Affero General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // This program is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU Affero General Public License for more details. // // You should have received a copy of the GNU Affero General Public License // along with this program. If not, see . // // All names, logos, and references to "Deltares" are registered trademarks of // 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.Data.Geometry; using Deltares.DamEngine.Data.Standard; namespace Deltares.DamEngine.Data.Geotechnics; /// /// Helper class for SoilProfile2D and SurfaceLine2 /// public static class SoilProfile2DSurfaceLineHelper { /// /// Determines if the surface line is inside the soil profile. /// /// The surface-line. /// The soil profile 2D /// public static bool IsSurfaceLineInsideSoilProfile2D(SurfaceLine2 surfaceLine, SoilProfile2D soilProfile2D) { const int precisionDecimals = 3; if (soilProfile2D == null) { return false; } if (Math.Round(surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelOutside).X, precisionDecimals) < Math.Round(soilProfile2D.Geometry.GetGeometryBounds().Left, precisionDecimals)) { return false; } if (Math.Round(surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelInside).X, precisionDecimals) > Math.Round(soilProfile2D.Geometry.GetGeometryBounds().Right, precisionDecimals)) { return false; } return true; } /// /// Combines the surface line with SoilProfile2D to create a new profile2D. /// /// The surface line. /// The SoilProfile2D. /// The default soil. /// The required shift of the geometry so it will fit the surface line /// public static SoilProfile2D CombineSurfaceLineWithSoilProfile2D(GeometryPointString surfaceLine, SoilProfile2D soilProfile2D, Soil defaultSoil, double shift) { if (soilProfile2D == null || soilProfile2D.Surfaces.Count == 0) return null; if (surfaceLine == null || surfaceLine.Points.Count == 0) return null; var result = new SoilProfile2D(); SoilProfile2D clonedProfile = soilProfile2D.Clone(); GeometryPointString clonedSurfaceLine = surfaceLine.Clone(); RoundCoordinates(clonedSurfaceLine, clonedProfile); shift = GeometryObject.RoundValue(shift); if (Math.Abs(shift) >= 0.0) { foreach (Point2D point in clonedProfile.Geometry.Points) point.X += shift; clonedProfile.Geometry.Rebox(); } if (clonedProfile.Geometry.Right < surfaceLine.GetMinX() || clonedProfile.Geometry.Left > surfaceLine.GetMaxX()) { throw new ArgumentException("Surface line is beyond the profile."); } // As the original profile may have been moved as cloned profile, a new clone is needed to perform all reset actions with. // The cloned profile is used to determine the original surfaces and preconsolidations. SoilProfile2D soilProfile2D2 = clonedProfile.Clone(); var oldSurfaces = new List(); oldSurfaces.AddRange((IEnumerable) soilProfile2D2.Surfaces); AddGeometryIntersectionsToSurfaceLine(clonedProfile.Geometry, ref clonedSurfaceLine); BuildGeometryModel(clonedSurfaceLine, clonedProfile, ref result); RemoveGeometryDataOfSoilProfileAboveSurfaceLine(clonedSurfaceLine, ref result); ReconstructSurfaces(ref result, clonedProfile, oldSurfaces, defaultSoil); ReconstructPreConsolidations(ref result, clonedProfile, shift); result.Geometry.Rebox(); result.Geometry.UpdateSurfaceLine(); return result; } /// /// Check if all the points of the surface line are above the bottom of the soil profile. /// /// /// /// public static bool IsSurfaceLineAboveBottomSoilProfile2D(SurfaceLine2 surfaceLine, SoilProfile2D soilProfile2D) { double bottom = soilProfile2D.Geometry.Bottom; foreach (GeometryPoint point in surfaceLine.Geometry.Points) { if (point.Z < bottom) { return false; } } return true; } private static void RoundCoordinates(GeometryPointString surfaceLine, SoilProfile2D soilProfile2D) { const double accuracy = 1e-7; foreach (Point2D point in soilProfile2D.Geometry.Points) { point.X = Point2D.RoundValue(point.X, accuracy); point.Z = Point2D.RoundValue(point.Z, accuracy); } soilProfile2D.Geometry.Right = GeometryObject.RoundValue(soilProfile2D.Geometry.Right, accuracy); soilProfile2D.Geometry.Left = GeometryObject.RoundValue(soilProfile2D.Geometry.Left, accuracy); soilProfile2D.Geometry.Bottom = GeometryObject.RoundValue(soilProfile2D.Geometry.Bottom, accuracy); surfaceLine.RoundPointsCoordinates(accuracy); } private static void AddGeometryIntersectionsToSurfaceLine(GeometryData geometry, ref GeometryPointString clonedSurfaceLine) { var pointList = new List(); foreach (GeometryLoop outerLoop in geometry.Surfaces.Select(surface => surface.OuterLoop)) { outerLoop.SyncCalcPoints(); pointList.AddRange(clonedSurfaceLine.IntersectionXzPointsWithGeometryPointList(outerLoop.CalcPoints)); } foreach (Point2D point in pointList) { if (clonedSurfaceLine.GetPointAt(point.X, point.Z) == null) clonedSurfaceLine.CalcPoints.Add(point); } clonedSurfaceLine.SyncPoints(); clonedSurfaceLine.SortPointsByXAscending(); } private static void BuildGeometryModel(GeometryPointString surfaceLine, SoilProfile2D soilProfile2D, ref SoilProfile2D result) { SoilProfile2D clonedProfile2 = soilProfile2D.Clone(); result.Geometry = clonedProfile2.Geometry; if (result.Geometry.Curves == null || result.Geometry.Curves.Count < 3) return; result.Geometry.Rebox(); if (surfaceLine == null || surfaceLine.Points.Count <= 1) return; AdaptGeometryToSurfaceLineLimits(surfaceLine, ref result); AddLeadingSurfaceLineToGeometry(surfaceLine, ref result); result.Geometry.NewlyEffectedPoints.AddRange((IEnumerable) result.Geometry.Points); result.Geometry.NewlyEffectedCurves.AddRange((IEnumerable) result.Geometry.Curves); result.Geometry.Surfaces.Clear(); result.Geometry.Curves.Clear(); result.Geometry.Points.Clear(); result.Geometry.RegenerateGeometry(); result.Geometry.DeleteLoosePoints(); result.Geometry.DeleteLooseCurves(); } private static void AdaptGeometryToSurfaceLineLimits(GeometryPointString surfaceLine, ref SoilProfile2D result) { GeometryBounds geometryBounds = surfaceLine.GetGeometryBounds(); if (result.Geometry.Left > geometryBounds.Left) GeometryHelper.ExtendGeometryLeft(result.Geometry, geometryBounds.Left); else if (result.Geometry.Left < geometryBounds.Left) GeometryHelper.CutGeometryLeft(result.Geometry, geometryBounds.Left); if (result.Geometry.Right > geometryBounds.Right) { GeometryHelper.CutGeometryRight(result.Geometry, geometryBounds.Right); } else { if (result.Geometry.Right < geometryBounds.Right) GeometryHelper.ExtendGeometryRight(result.Geometry, geometryBounds.Right); } } private static void AddLeadingSurfaceLineToGeometry(GeometryPointString surfaceLine, ref SoilProfile2D result) { surfaceLine.SyncCalcPoints(); IList points = surfaceLine.CalcPoints; GeometryData geometry = result.Geometry; Point2D current1 = points[0].Clone(); var existingPoint1 = geometry.GetPointAtLocation(current1, GeometryConstants.Accuracy); if (existingPoint1 == null) { Point2D[] leftPoints = geometry.GetLeftPoints().OrderBy(x => x.Z).ToArray(); geometry.NewlyEffectedPoints.Add(current1); if (leftPoints[^1].Z < current1.Z) { GeometryCurve newLeftCurve = new GeometryCurve(leftPoints[^1], current1); result.Geometry.NewlyEffectedCurves.Add(newLeftCurve); } } else { current1 = existingPoint1; } Point2D current2; int index = 1; while (index < points.Count) { current2 = points[index].Clone(); var existingPoint2 = geometry.GetPointAtLocation(current2, GeometryConstants.Accuracy); if (existingPoint2 == null) { if (index == checked (points.Count - 1)) { Point2D[] rightPoints = geometry.GetRightPoints().OrderBy(x => x.Z).ToArray(); geometry.NewlyEffectedPoints.Add(current2); if (rightPoints[^1].Z < current2.Z) { GeometryCurve newRightCurve = new GeometryCurve(rightPoints[^1], current2); result.Geometry.NewlyEffectedCurves.Add(newRightCurve); } } else { geometry.NewlyEffectedPoints.Add(current2); } } else { current2 = existingPoint2; } GeometryCurve newCurve = new GeometryCurve(current1, current2); if (GeometryHelper.DoesCurveExist(result.Geometry, newCurve) == null) { result.Geometry.NewlyEffectedCurves.Add(newCurve); } current1 = current2; index++; } } private static void RemoveGeometryDataOfSoilProfileAboveSurfaceLine(GeometryPointString surfaceLine, ref SoilProfile2D result) { GeometryData geometry = result.Geometry; RemoveGeometryDataAboveSurfaceLine(surfaceLine, ref geometry); result.Geometry = geometry; } private static void RemoveGeometryDataAboveSurfaceLine(GeometryPointString surfaceLine, ref GeometryData result) { if (surfaceLine == null || surfaceLine.Points.Count <= 1) return; foreach (GeometrySurface geometrySurface in result.Surfaces.ToArray()) { if (surfaceLine.GetPosition(geometrySurface.OuterLoop) == RelativeXzPosition.AboveGeometricLine) { result.Surfaces.Remove(geometrySurface); foreach (GeometryCurve curve in geometrySurface.OuterLoop.CurveList) { if (curve.SurfaceAtLeft == geometrySurface) curve.SurfaceAtLeft = null; if (curve.SurfaceAtRight == geometrySurface) curve.SurfaceAtRight = null; } foreach (GeometryLoop geometryLoop in result.Loops.ToArray()) { GeometryLoop loop = geometryLoop; if (!result.Surfaces.Any((Func)(s => s.OuterLoop == loop)) && !result.Surfaces.Any((Func)(s => s.InnerLoops.Contains(loop)))) { result.Loops.Remove(loop); } } } } result.DeleteLooseCurves(); result.DeleteLoosePoints(); } private static void ReconstructSurfaces(ref SoilProfile2D result, SoilProfile2D soilProfile2D, List oldSurfaces, Soil defaultSoil) { result.Surfaces.Clear(); result.Name = soilProfile2D.Name; foreach (GeometrySurface surface in result.Geometry.Surfaces) { var soilLayer2D1 = new SoilLayer2D { GeometrySurface = surface, Soil = defaultSoil }; SoilLayer2D soilLayer2D2 = soilLayer2D1; SoilLayer2D layerFromOldSurfaces = SoilProfile2D.GetOriginalLayerFromOldSurfaces(surface, (IEnumerable) oldSurfaces); if (layerFromOldSurfaces != null && layerFromOldSurfaces.Soil != null) { soilLayer2D2.Soil = layerFromOldSurfaces.Soil; soilLayer2D2.IsAquifer = layerFromOldSurfaces.IsAquifer; soilLayer2D2.WaterpressureInterpolationModel = layerFromOldSurfaces.WaterpressureInterpolationModel; } if (layerFromOldSurfaces == null) { SoilLayer2D oldLayer = null; bool isDefaultLayer = IsLayerAboveOriginalSurfaceLine(soilLayer2D2, soilProfile2D.Geometry.SurfaceLine); if (!isDefaultLayer) { oldLayer = GetLayerIfSurfaceIsLeftOrRightOfOldSurfaces(surface, oldSurfaces); isDefaultLayer = oldLayer == null; } if (!isDefaultLayer) { soilLayer2D2.IsAquifer = oldLayer.IsAquifer; soilLayer2D2.WaterpressureInterpolationModel = oldLayer.WaterpressureInterpolationModel; soilLayer2D2.SoilName = oldLayer.SoilName; soilLayer2D2.Soil = oldLayer.Soil; } else { soilLayer2D2.IsAquifer = false; soilLayer2D2.WaterpressureInterpolationModel = WaterpressureInterpolationModel.Automatic; soilLayer2D2.SoilName = defaultSoil.Name; soilLayer2D2.Soil = defaultSoil; } } result.Surfaces.Add(soilLayer2D2); } } internal static bool IsLayerAboveOriginalSurfaceLine(SoilLayer2D layer, GeometryPointString surfaceLine) { if (layer == null || surfaceLine == null) { return false; } return layer.GeometrySurface.OuterLoop.Points.All(point => !point.Z.IsLessThan(surfaceLine.GetZatX(point.X))); } private static void ReconstructPreConsolidations(ref SoilProfile2D result, SoilProfile2D cloneProfile, double shift) { foreach (PreConsolidationStress preconsolidationStress in cloneProfile.PreconsolidationStresses) { preconsolidationStress.X += shift; result.PreconsolidationStresses.Add(preconsolidationStress); } } private static SoilLayer2D GetLayerIfSurfaceIsLeftOrRightOfOldSurfaces(GeometrySurface surface, List oldSurfaces) { double[] xMin = surface.OuterLoop.CalcPoints.Select(p => p.X).OrderBy(x => x).Distinct().ToArray(); double[] zMin = surface.OuterLoop.CalcPoints.Select(p => p.Z).OrderBy(z => z).Distinct().ToArray(); var leftPoint = new Point2D { X = xMin[0] - 0.1, Z = zMin[0] + 0.1 }; SoilLayer2D soilLayer = FindLayerFromPointInOldSurfaces(oldSurfaces, leftPoint); if (soilLayer != null) { return soilLayer; } var rightPoint = new Point2D { X = xMin[^1] + 0.1, Z = zMin[0] + 0.1 }; soilLayer = FindLayerFromPointInOldSurfaces(oldSurfaces, rightPoint); return soilLayer; } private static SoilLayer2D FindLayerFromPointInOldSurfaces(List oldSurfaces, Point2D searchPoint) { for (var i = 0; i < oldSurfaces.Count; i++) { bool find = oldSurfaces[i].GeometrySurface.OuterLoop.IsPointInLoopArea(searchPoint); if (find) { return oldSurfaces[i]; } } return null; } }