Index: DamEngine/trunk/src/Deltares.DamEngine.Data.Tests/Geotechnics/SoilSurfaceProfileTests.cs =================================================================== diff -u -r3247 -r3254 --- DamEngine/trunk/src/Deltares.DamEngine.Data.Tests/Geotechnics/SoilSurfaceProfileTests.cs (.../SoilSurfaceProfileTests.cs) (revision 3247) +++ DamEngine/trunk/src/Deltares.DamEngine.Data.Tests/Geotechnics/SoilSurfaceProfileTests.cs (.../SoilSurfaceProfileTests.cs) (revision 3254) @@ -328,6 +328,197 @@ }, bottomSoilLayer2D.GeometrySurface.OuterLoop.CurveList); } + [Test] + public void ConvertToSoilProfile2D_WithSurfaceLineIntersectingSoilLayers_ReturnsExpectedSoilProfile() + { + // Setup + const string bottomLayerName = "BottomLayer"; + const string middleLayerName = "MiddleLayer"; + const string topLayerName = "TopLayer"; + + SoilLayer1D bottomLayer = CreateSoilLayer(-1.25, bottomLayerName); + SoilLayer1D middleLayer = CreateSoilLayer(0, middleLayerName); + + var profile = new SoilProfile1D(); + profile.Layers.Add(bottomLayer); + profile.Layers.Add(middleLayer); + profile.BottomLevel = -10; + + SurfaceLine2 surfaceLine = CreateSurfaceLine(new[] + { + new GeometryPoint(0, 2.5), + new GeometryPoint(5, -2.5), + new GeometryPoint(10, 2.5) + }); + + var soilSurfaceProfile = new SoilSurfaceProfile + { + SoilProfile = profile, + SurfaceLine2 = surfaceLine, + DikeEmbankmentMaterial = new Soil + { + Name = topLayerName + } + }; + + // Call + SoilProfile2D soilProfile2D = soilSurfaceProfile.ConvertToSoilProfile2D(); + + // Assert + var soilLayer2Ds = soilProfile2D.Surfaces; + Assert.That(soilLayer2Ds, Has.Count.EqualTo(5)); + + var topSoilLayerSurfaces = soilLayer2Ds.Where(l => string.Equals(l.Name, topLayerName)).ToArray(); + Assert.That(topSoilLayerSurfaces.Length, Is.EqualTo(2)); + Assert.That(topSoilLayerSurfaces.All(s => ReferenceEquals(s.Soil, soilSurfaceProfile.DikeEmbankmentMaterial)), Is.True); + AssertGeometry(new [] + { + new GeometryCurve(new Point2D(0, 2.5), new Point2D(2.5, 0)), + new GeometryCurve(new Point2D(2.5, 0), new Point2D(0, 0)), + new GeometryCurve(new Point2D(0, 0), new Point2D(0, 2.5)) + }, topSoilLayerSurfaces[0].GeometrySurface.OuterLoop.CurveList); + AssertGeometry(new[] + { + new GeometryCurve(new Point2D(7.5, 0), new Point2D(10, 2.5)), + new GeometryCurve(new Point2D(10, 2.5), new Point2D(10, 0)), + new GeometryCurve(new Point2D(10, 0), new Point2D(7.5, 0)) + }, topSoilLayerSurfaces[1].GeometrySurface.OuterLoop.CurveList); + + var middleLayerSurfaces = soilLayer2Ds.Where(l => string.Equals(l.Name, middleLayerName)).ToArray(); + Assert.That(middleLayerSurfaces.Length, Is.EqualTo(2)); + Assert.That(middleLayerSurfaces.All(s => ReferenceEquals(s.Soil, middleLayer.Soil)), Is.True); + AssertGeometry(new[] + { + new GeometryCurve(new Point2D(0, 0), new Point2D(2.5, 0)), + new GeometryCurve(new Point2D(2.5, 0), new Point2D(3.75, -1.25)), + new GeometryCurve(new Point2D(3.75, -1.25), new Point2D(2.5, -1.25)), + new GeometryCurve(new Point2D(2.5, -1.25), new Point2D(0, -1.25)), + new GeometryCurve(new Point2D(0, -1.25), new Point2D(0, 0)) + }, middleLayerSurfaces[0].GeometrySurface.OuterLoop.CurveList); + AssertGeometry(new[] + { + new GeometryCurve(new Point2D(7.5, 0), new Point2D(10, 0)), + new GeometryCurve(new Point2D(10, 0), new Point2D(10, -1.25)), + new GeometryCurve(new Point2D(10, -1.25), new Point2D(7.5, -1.25)), + new GeometryCurve(new Point2D(7.5, -1.25), new Point2D(6.25, -1.25)), + new GeometryCurve(new Point2D(6.25, -1.25), new Point2D(7.5, 0)) + }, middleLayerSurfaces[1].GeometrySurface.OuterLoop.CurveList); + + SoilLayer2D bottomLayerSurface = soilLayer2Ds.Single(l => string.Equals(l.Name, bottomLayerName)); + Assert.That(bottomLayerSurface.Soil, Is.SameAs(bottomLayer.Soil)); + AssertGeometry(new[] + { + new GeometryCurve(new Point2D(0, -1.25), new Point2D(2.5, -1.25)), + new GeometryCurve(new Point2D(2.5, -1.25), new Point2D(3.75, -1.25)), + new GeometryCurve(new Point2D(3.75, -1.25), new Point2D(5, -2.5)), + new GeometryCurve(new Point2D(5, -2.5), new Point2D(6.25, -1.25)), + new GeometryCurve(new Point2D(6.25, -1.25), new Point2D(7.5, -1.25)), + new GeometryCurve(new Point2D(7.5, -1.25), new Point2D(10, -1.25)), + new GeometryCurve(new Point2D(10, -1.25), new Point2D(10, -10)), + new GeometryCurve(new Point2D(10, -10), new Point2D(7.5, -10)), + new GeometryCurve(new Point2D(7.5, -10), new Point2D(6.25, -10)), + new GeometryCurve(new Point2D(6.25, -10), new Point2D(5, -10)), + new GeometryCurve(new Point2D(5, -10), new Point2D(3.75, -10)), + new GeometryCurve(new Point2D(3.75, -10), new Point2D(2.5, -10)), + new GeometryCurve(new Point2D(2.5, -10), new Point2D(0, -10)), + new GeometryCurve(new Point2D(0, -10), new Point2D(0, -1.25)) + + }, bottomLayerSurface.GeometrySurface.OuterLoop.CurveList); + } + + [Test] + public void ConvertToSoilProfile2D_WithSurfaceLineIntersectingSoilLayerAndInflectionAtBottomIntersection_ReturnsExpectedSoilProfile() + { + // Setup + const string bottomLayerName = "BottomLayer"; + const string middleLayerName = "MiddleLayer"; + const string topLayerName = "TopLayer"; + + SoilLayer1D bottomLayer = CreateSoilLayer(-2.5, bottomLayerName); + SoilLayer1D middleLayer = CreateSoilLayer(0, middleLayerName); + + var profile = new SoilProfile1D(); + profile.Layers.Add(bottomLayer); + profile.Layers.Add(middleLayer); + profile.BottomLevel = -10; + + SurfaceLine2 surfaceLine = CreateSurfaceLine(new[] + { + new GeometryPoint(0, 2.5), + new GeometryPoint(5, -2.5), + new GeometryPoint(10, 2.5) + }); + + var soilSurfaceProfile = new SoilSurfaceProfile + { + SoilProfile = profile, + SurfaceLine2 = surfaceLine, + DikeEmbankmentMaterial = new Soil + { + Name = topLayerName + } + }; + + // Call + SoilProfile2D soilProfile2D = soilSurfaceProfile.ConvertToSoilProfile2D(); + + // Assert + var soilLayer2Ds = soilProfile2D.Surfaces; + Assert.That(soilLayer2Ds, Has.Count.EqualTo(5)); + + var topSoilLayerSurfaces = soilLayer2Ds.Where(l => string.Equals(l.Name, topLayerName)).ToArray(); + Assert.That(topSoilLayerSurfaces.Length, Is.EqualTo(2)); + Assert.That(topSoilLayerSurfaces.All(s => ReferenceEquals(s.Soil, soilSurfaceProfile.DikeEmbankmentMaterial)), Is.True); + AssertGeometry(new[] + { + new GeometryCurve(new Point2D(0, 2.5), new Point2D(2.5, 0)), + new GeometryCurve(new Point2D(2.5, 0), new Point2D(0, 0)), + new GeometryCurve(new Point2D(0, 0), new Point2D(0, 2.5)) + }, topSoilLayerSurfaces[0].GeometrySurface.OuterLoop.CurveList); + AssertGeometry(new[] + { + new GeometryCurve(new Point2D(7.5, 0), new Point2D(10, 2.5)), + new GeometryCurve(new Point2D(10, 2.5), new Point2D(10, 0)), + new GeometryCurve(new Point2D(10, 0), new Point2D(7.5, 0)) + }, topSoilLayerSurfaces[1].GeometrySurface.OuterLoop.CurveList); + + var middleLayerSurfaces = soilLayer2Ds.Where(l => string.Equals(l.Name, middleLayerName)).ToArray(); + Assert.That(middleLayerSurfaces.Length, Is.EqualTo(2)); + Assert.That(middleLayerSurfaces.All(s => ReferenceEquals(s.Soil, middleLayer.Soil)), Is.True); + AssertGeometry(new[] + { + new GeometryCurve(new Point2D(0, 0), new Point2D(2.5, 0)), + new GeometryCurve(new Point2D(2.5, 0), new Point2D(5, -2.5)), + new GeometryCurve(new Point2D(5, -2.5), new Point2D(2.5, -2.5)), + new GeometryCurve(new Point2D(2.5, -2.5), new Point2D(0, -2.5)), + new GeometryCurve(new Point2D(0, -2.5), new Point2D(0, 0)) + }, middleLayerSurfaces[0].GeometrySurface.OuterLoop.CurveList); + AssertGeometry(new[] + { + new GeometryCurve(new Point2D(7.5, 0), new Point2D(10, 0)), + new GeometryCurve(new Point2D(10, 0), new Point2D(10, -2.5)), + new GeometryCurve(new Point2D(10, -2.5), new Point2D(7.5, -2.5)), + new GeometryCurve(new Point2D(7.5, -2.5), new Point2D(5, -2.5)), + new GeometryCurve(new Point2D(5, -2.5), new Point2D(7.5, 0)) + }, middleLayerSurfaces[1].GeometrySurface.OuterLoop.CurveList); + + SoilLayer2D bottomLayerSurface = soilLayer2Ds.Single(l => string.Equals(l.Name, bottomLayerName)); + Assert.That(bottomLayerSurface.Soil, Is.SameAs(bottomLayer.Soil)); + AssertGeometry(new[] + { + new GeometryCurve(new Point2D(0, -2.5), new Point2D(2.5, -2.5)), + new GeometryCurve(new Point2D(2.5, -2.5), new Point2D(5, -2.5)), + new GeometryCurve(new Point2D(5, -2.5), new Point2D(7.5, -2.5)), + new GeometryCurve(new Point2D(7.5, -2.5), new Point2D(10, -2.5)), + new GeometryCurve(new Point2D(10, -2.5), new Point2D(10, -10)), + new GeometryCurve(new Point2D(10, -10), new Point2D(7.5, -10)), + new GeometryCurve(new Point2D(7.5, -10), new Point2D(5, -10)), + new GeometryCurve(new Point2D(5, -10), new Point2D(2.5, -10)), + new GeometryCurve(new Point2D(2.5, -10), new Point2D(0, -10)), + new GeometryCurve(new Point2D(0, -10), new Point2D(0, -2.5)) + }, bottomLayerSurface.GeometrySurface.OuterLoop.CurveList); + } + private static void AssertGeometry(IEnumerable expectedCurves, IEnumerable actualCurves) { int nrOfExpectedCurves = expectedCurves.Count(); Index: DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/SoilSurfaceProfile.cs =================================================================== diff -u -r3248 -r3254 --- DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/SoilSurfaceProfile.cs (.../SoilSurfaceProfile.cs) (revision 3248) +++ DamEngine/trunk/src/Deltares.DamEngine.Data/Geotechnics/SoilSurfaceProfile.cs (.../SoilSurfaceProfile.cs) (revision 3254) @@ -437,15 +437,27 @@ private static IEnumerable CreateSoilLayerWithSurfaceLineIntersection(SoilLayer1D soilLayer, SurfaceLine2 surfaceLine, IEnumerable xCoordinates) { + // Surface cannot be determined if there are not at least 2 coordinates + var xCoordinateArray = xCoordinates.ToArray(); + if (xCoordinateArray.Length < 2) + { + return Enumerable.Empty(); + } + double layerTopLevel = soilLayer.TopLevel; double layerBottomLevel = soilLayer.BottomLevel; // Determine the points that are part of the geometry var currentArea = new EnclosedArea(); var enclosedAreas = new List(); + double previousZCoordinate = double.NegativeInfinity; - foreach (double xCoordinate in xCoordinates) + int j = 1; + for (int i = 0; i < xCoordinateArray.Length; i++) { + double xCoordinate = xCoordinateArray[i]; + + double nextZCoordinate = surfaceLine.Geometry.GetZatX(xCoordinateArray[j]); double zCoordinateSurfaceLine = surfaceLine.Geometry.GetZatX(xCoordinate); double currentZCoordinate = zCoordinateSurfaceLine; if (zCoordinateSurfaceLine > layerTopLevel) @@ -459,15 +471,27 @@ } // Determine if a new area should be started. This should only happen when: - if (previousZCoordinate > layerBottomLevel // The surface line already starts above the soil layer bottom - && currentZCoordinate < previousZCoordinate // The surface line has a decreasing trend - && currentZCoordinate < layerBottomLevel) // The surface line crossed the bottom of the layer + if (previousZCoordinate > layerBottomLevel // The area already started above or at the soil layer bottom + && zCoordinateSurfaceLine < previousZCoordinate // The surface line has a decreasing trend wrt the previous z coordinate + && zCoordinateSurfaceLine <= layerBottomLevel) // The surface line crossed the bottom of the layer { enclosedAreas.Add(currentArea); currentArea = new EnclosedArea(); + + // If the next z coordinate is increasing, it means that the surface line is intersecting + // the bottom level with an inflection point. In this situation, the current coordinate should be + // added to the new area + if (nextZCoordinate > currentZCoordinate) + { + currentArea.AddCoordinate(new Point2D(xCoordinate, currentZCoordinate)); + } + } previousZCoordinate = currentZCoordinate; + + // End of array, take the last Z coordinate as the next coordinate + j = i == xCoordinateArray.Length - 2 ? j : j + 1; } // The last area needs to be added manually @@ -506,7 +530,7 @@ } // Close the geometry by connecting the first and last points with each other - // This closes the outerloop in a clockwise definition. + // This closes the outer loop in a clockwise definition. curves.Add(new GeometryCurve(pointsToConvert.Last(), pointsToConvert.First())); return new LayerData(pointsToConvert, curves, soilLayer.Soil, soilLayer.IsAquifer, soilLayer.WaterpressureInterpolationModel); @@ -624,8 +648,6 @@ var xCoordinatesToTraverse = GetXCoordinates(filteredSoilLayers, surfaceLine2).ToArray(); var layerData = CreateLayerData(filteredSoilLayers, surfaceLine2, xCoordinatesToTraverse); - - // Build geometry BuildGeometryModel(Geometry, layerData, originalSoilProfile1D, surfaceLine2); BuildSoilLayer2D(Surfaces, layerData, Geometry);