// 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.GetPoint2D(CharacteristicPointType.SurfaceLevelOutside).X, precisionDecimals) < Math.Round(soilProfile2D.Geometry.GetGeometryBounds().Left, precisionDecimals)) { return false; } if (Math.Round(surfaceLine.CharacteristicPoints.GetPoint2D(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 dike embankment soil. /// public static SoilProfile2D CombineSurfaceLineWithSoilProfile2D(GeometryPointString surfaceLine, SoilProfile2D soilProfile2D, Soil defaultSoil) { if (soilProfile2D == null || soilProfile2D.Surfaces.Count == 0) return null; if (surfaceLine == null || surfaceLine.Points.Count == 0) return null; SoilProfile2D clonedProfile = soilProfile2D.Clone(); GeometryPointString clonedSurfaceLine = surfaceLine.Clone(); RoundCoordinates(clonedProfile.Geometry, clonedSurfaceLine); if (clonedProfile.Geometry.Right < surfaceLine.GetMinX() || clonedProfile.Geometry.Left > surfaceLine.GetMaxX()) { throw new ArgumentException("Surface line is beyond the profile."); } var oldSurfaces = new List(); oldSurfaces.AddRange((IEnumerable) soilProfile2D.Surfaces); var result = new SoilProfile2D(); result.Name = soilProfile2D.Name; result.Geometry = CreateNewGeometryForSoilProfile2DByCombiningItsGeometryWithSurfaceLine(clonedSurfaceLine, clonedProfile.Geometry); RemoveGeometryDataOfSoilProfileAboveSurfaceLine(clonedSurfaceLine, ref result); ReconstructSurfaces(ref result, clonedSurfaceLine, oldSurfaces, defaultSoil); ReconstructPreConsolidations(ref result, clonedProfile); result.Geometry.Rebox(); result.Geometry.UpdateSurfaceLine(); RoundCoordinates(result.Geometry); 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 (Point2D point in surfaceLine.Geometry.Points) { if (point.Z < bottom) { return false; } } return true; } private static void RoundCoordinates(GeometryData geometry, GeometryPointString surfaceLine = null) { foreach (Point2D point in geometry.Points) { point.X = Point2D.RoundValue(point.X); point.Z = Point2D.RoundValue(point.Z); } geometry.Right = GeometryObject.RoundValue(geometry.Right); geometry.Left = GeometryObject.RoundValue(geometry.Left); geometry.Bottom = GeometryObject.RoundValue(geometry.Bottom); if (surfaceLine != null) { surfaceLine.RoundPointsCoordinates(GeometryConstants.Accuracy); } } private static GeometryData CreateNewGeometryForSoilProfile2DByCombiningItsGeometryWithSurfaceLine( GeometryPointString surfaceLine, GeometryData originalGeometry) { var resultGeometry = originalGeometry.Clone(); if (resultGeometry == null) { return null; } if (resultGeometry.Curves == null || resultGeometry.Curves.Count < 3) { return null; } resultGeometry.Rebox(); if (surfaceLine == null || surfaceLine.Points.Count <= 2) { return null; } if (resultGeometry.NewlyEffectedCurves.Count > 0 || resultGeometry.NewlyEffectedPoints.Count > 0) { return null; } AdaptGeometryToSurfaceLineLimits(surfaceLine, ref resultGeometry); AddLeadingSurfaceLineToGeometry(surfaceLine, ref resultGeometry); resultGeometry.NewlyEffectedPoints.AddRange((IEnumerable) resultGeometry.Points); resultGeometry.NewlyEffectedCurves.AddRange((IEnumerable) resultGeometry.Curves); resultGeometry.Surfaces.Clear(); resultGeometry.Curves.Clear(); resultGeometry.Points.Clear(); resultGeometry.RegenerateGeometry(); resultGeometry.DeleteLoosePoints(); resultGeometry.DeleteLooseCurves(); return resultGeometry; } private static void AdaptGeometryToSurfaceLineLimits(GeometryPointString surfaceLine, ref GeometryData geometryData) { GeometryBounds geometryBounds = surfaceLine.GetGeometryBounds(); if (geometryBounds.Left.IsLessThan(geometryData.Left)) GeometryHelper.ExtendGeometryLeft(geometryData, geometryBounds.Left); else if (geometryData.Left.IsLessThan(geometryBounds.Left)) GeometryHelper.CutGeometryLeft(geometryData, geometryBounds.Left); if (geometryBounds.Right.IsLessThan(geometryData.Right)) { GeometryHelper.CutGeometryRight(geometryData, geometryBounds.Right); } else { if (geometryData.Right.IsLessThan(geometryBounds.Right)) GeometryHelper.ExtendGeometryRight(geometryData, geometryBounds.Right); } } private static void AddLeadingSurfaceLineToGeometry(GeometryPointString surfaceLine, ref GeometryData geometryData) { geometryData.NewlyEffectedCurves.Clear(); geometryData.NewlyEffectedPoints.Clear(); IList points = surfaceLine.Points; GeometryData geometry = geometryData; Point2D currentPoint1 = points[0].Clone(); var existingPoint1 = geometry.GetPointAtLocation(currentPoint1); if (existingPoint1 == null) { Point2D[] leftPoints = geometry.GetLeftPoints().OrderBy(x => x.Z).ToArray(); geometry.NewlyEffectedPoints.Add(currentPoint1); if (leftPoints[^1].Z < currentPoint1.Z) { // When the first point of the surface line is above the geometry (left points) add a curve to the geometry // from the last (= top) point of the left points to the first point of the surface line. GeometryCurve newLeftCurve = new GeometryCurve(leftPoints[^1], currentPoint1); geometryData.NewlyEffectedCurves.Add(newLeftCurve); } } else { currentPoint1 = existingPoint1; } Point2D currentPoint2; int index = 1; while (index < points.Count) { currentPoint2 = points[index].Clone(); var existingPoint2 = geometry.GetPointAtLocation(currentPoint2); if (existingPoint2 == null) { if (index == checked (points.Count - 1)) { Point2D[] rightPoints = geometry.GetRightPoints().OrderBy(x => x.Z).ToArray(); geometry.NewlyEffectedPoints.Add(currentPoint2); if (rightPoints[^1].Z < currentPoint2.Z) { // When the last point of the surface line is above the geometry (right points) add a curve to the // geometry from the last (= top) point of the right points to the last point of the surface line. GeometryCurve newRightCurve = new GeometryCurve(rightPoints[^1], currentPoint2); geometryData.NewlyEffectedCurves.Add(newRightCurve); } } else { geometry.NewlyEffectedPoints.Add(currentPoint2); } } else { currentPoint2 = existingPoint2; } GeometryCurve newCurve = new GeometryCurve(currentPoint1, currentPoint2); if (GeometryHelper.DoesCurveExist(geometryData, newCurve) == null) { geometryData.NewlyEffectedCurves.Add(newCurve); } currentPoint1 = currentPoint2; 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, GeometryPointString surfaceLine, List oldSurfaces, Soil defaultSoil) { result.Surfaces.Clear(); foreach (GeometrySurface surface in result.Geometry.Surfaces) { var soilLayer2D1 = new SoilLayer2D { GeometrySurface = surface, Soil = defaultSoil }; SoilLayer2D soilLayer2D2 = soilLayer2D1; SoilLayer2D layerFromOldSurfaces = SoilProfile2D.DetermineOriginalLayerFromOldSurfaces(surface, 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, surfaceLine); if (!isDefaultLayer) { oldLayer = DetermineLayerIfSurfaceIsLeftOrRightOfOldSurfaces(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 clonedProfile) { foreach (PreConsolidationStress preconsolidationStress in clonedProfile.PreconsolidationStresses) { result.PreconsolidationStresses.Add(preconsolidationStress); } } private static SoilLayer2D DetermineLayerIfSurfaceIsLeftOrRightOfOldSurfaces(GeometrySurface surface, List oldSurfaces) { double[] xMin = surface.OuterLoop.Points.Select(p => p.X).OrderBy(x => x).Distinct().ToArray(); double[] zMin = surface.OuterLoop.Points.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; } }