//-----------------------------------------------------------------------
//
// Copyright (c) 2009 Deltares. All rights reserved.
//
// B.S.T.I.M. The
// tom.the@deltares.nl
// 18-05-2009
// Contains class to create PL-Line
//-----------------------------------------------------------------------
using System.Runtime.Serialization;
using Deltares.Geometry;
using Deltares.Geotechnics.Soils;
using Deltares.Geotechnics.SurfaceLines;
using Deltares.Mathematics;
using Deltares.Standard.Extensions;
using Deltares.Uplift;
namespace Deltares.Dam.Data
{
using System;
using System.Collections.Generic;
using System.Linq;
using Deltares.Geotechnics;
using Deltares.Soilbase;
[Serializable]
public class PLLinesCreatorException : Exception
{
public PLLinesCreatorException()
{
}
public PLLinesCreatorException(string message)
: base(message)
{
}
public PLLinesCreatorException(string message, Exception inner)
: base(message, inner)
{
}
protected PLLinesCreatorException(SerializationInfo info, StreamingContext context)
: base(info, context)
{
}
}
///
///
///
public class PLLinesCreator
{
private const double cUpliftFactorEquilibrium = 1.0;
private const double cOffsetPhreaticLineBelowSurface = 0.01;
private ModelParametersForPLLines modelParametersForPLLines = new ModelParametersForPLLines();
public bool IsUseOvenDryUnitWeight { get; set; }
public SoilGeometryType SoilGeometryType { get; set; }
public SoilProfile1D SoilProfile { get; set; }
public string SoilGeometry2DName { get; set; }
public Soil DikeEmbankmentMaterial { get; set; }
public SoilbaseDB SoilBaseDB { get; set; }
public SoilList SoilList { get; set; }
public double? HeadInPLLine2 { get; set; }
private double? headInPLLine3 { get; set; }
private double? headInPLLine4 { get; set; }
public double XSoilGeometry2DOrigin { get; set; }
private bool isAdjustPL3AndPL4SoNoUpliftWillOccurEnabled = true;
private bool isHydraulicShortcut = false;
private double waterLevelPolder;
private double waterLevelRiverHigh;
private double? waterLevelRiverLow;
private bool isUseLowWaterLevel;
private SurfaceLine2 surfaceLine;
public bool IsAdjustPL3AndPL4SoNoUpliftWillOccurEnabled { get { return isAdjustPL3AndPL4SoNoUpliftWillOccurEnabled; } set { isAdjustPL3AndPL4SoNoUpliftWillOccurEnabled = value; } }
public double WaterLevelRiverHigh { get { return waterLevelRiverHigh; } set { waterLevelRiverHigh = value; } }
public double? WaterLevelRiverLow { get { return waterLevelRiverLow; } set { waterLevelRiverLow = value; } }
public bool IsUseLowWaterLevel { get { return isUseLowWaterLevel; } set { isUseLowWaterLevel = value; } }
public double WaterLevelPolder { get { return waterLevelPolder; } set { waterLevelPolder = value; } }
public double? HeadInPLLine3 { get { return headInPLLine3; } set { headInPLLine3 = value; } }
public double? HeadInPLLine4 { get { return headInPLLine4; } set { headInPLLine4 = value; } }
/// Aggregation relationship
public SurfaceLine2 SurfaceLine
{
get { return surfaceLine; }
set { surfaceLine = value; }
}
public IList GaugePLLines { get; set; }
public IList Gauges { get; set; }
public double GaugeMissVal { get; set; }
public double PlLineOffsetBelowDikeTopAtRiver { get; set; }
public double PlLineOffsetBelowDikeTopAtPolder { get; set; }
public virtual double PlLineOffsetBelowShoulderBaseInside { get; set; }
public virtual double PlLineOffsetBelowDikeToeAtPolder { get; set; }
public double? PlLineOffsetBelowDikeCrestMiddle { get; set; }
public double? PlLineOffsetFactorBelowShoulderCrest { get; set; }
public bool? UsePlLineOffsetBelowDikeCrestMiddle { get; set; }
public bool? UsePlLineOffsetFactorBelowShoulderCrest { get; set; }
public PhreaticAdaptionType? NWOPhreaticAdaption { get; set; }
private double WaterLevelToUse()
{
if (isUseLowWaterLevel)
{
return this.waterLevelRiverLow.Value;
}
else
{
return this.WaterLevelRiverHigh;
}
}
private double HeadPL3TakingInAccountHydraulicShortcut
{
get
{
double waterLevel = WaterLevelToUse();
// If no value known for headPl3 then use waterlevel for headPL3
if (HeadInPLLine3 == null)
{
return waterLevel;
}
// If hydraulic shortcut and no inbetween aquifer layer then taka into account hydraulic shortcut
// If hydraulic shortcut then use waterlevel for headPL3
if (IsHydraulicShortcut)
{
if (SoilGeometryType == SoilGeometryType.SoilGeometry1D)
{
if (SoilProfile.InBetweenAquiferLayer == null)
{
return waterLevel;
}
}
return HeadInPLLine3.Value;
}
else
{
return HeadInPLLine3.Value;
}
}
}
private double GetHeadPL4TakingInAccountHydraulicShortcut(DamType damType)
{
// If hydraulic shortcut or no value known for headPl4 then use waterlevel for headPL4
double waterLevel = WaterLevelToUse();
if (IsHydraulicShortcut)
{
// If hydraulic shortcut then use waterlevel for headPL4
return waterLevel;
}
else
{
if (HeadInPLLine4 == null)
{
// If no hydraulic shortcut and no value known for headPl4 then use polderlevel for headPL4
if (damType == DamType.Primary)
{
return waterLevel;
}
else
{
return waterLevelPolder;
}
}
else
{
// If no hydraulic shortcut and value is specified for headPl4 then use this value for headPL4
return HeadInPLLine4.Value;
}
}
}
private PLLine currentPL1Line = null; // is needed when calculating uplift reduction for PL3 and Pl4
// Output
private double pl3MinUplift;
private double pl3HeadAdjusted;
private double pl3LocationXMinUplift;
public double Pl3MinUplift { get { return pl3MinUplift; } }
public double Pl3HeadAdjusted { get { return pl3HeadAdjusted; } }
public double Pl3LocationXMinUplift { get { return pl3LocationXMinUplift; } }
private double pl4MinUplift;
private double pl4HeadAdjusted;
private double pl4LocationXMinUplift;
public double Pl4MinUplift { get { return pl4MinUplift; } }
public double Pl4HeadAdjusted { get { return pl4HeadAdjusted; } }
public double Pl4LocationXMinUplift { get { return pl4LocationXMinUplift; } }
///
/// Constructor
///
public PLLinesCreator()
{
SoilGeometryType = SoilGeometryType.SoilGeometry1D;
PlLineOffsetBelowDikeTopAtRiver = 0.5; // Default value
PlLineOffsetBelowDikeTopAtPolder = 1.5; // Default value
PlLineOffsetBelowShoulderBaseInside = 0.1; // Default value
PlLineOffsetBelowDikeToeAtPolder = 0.1; // Default value
IsUseOvenDryUnitWeight = false;
}
///
///
///
///
///
private SoilProfile1D GetSoilProfileBelowPoint(double xCoordinate)
{
switch (this.SoilGeometryType)
{
case SoilGeometryType.SoilGeometry1D:
var soilProfile = new SoilProfile1D();
soilProfile.Assign(this.SoilProfile);
return soilProfile;
case SoilGeometryType.SoilGeometry2D:
var geometry2DTo1DConverter = new Geometry2DTo1DConverter(this.SoilGeometry2DName, this.SurfaceLine, this.DikeEmbankmentMaterial, this.SoilBaseDB, this.SoilList, -this.XSoilGeometry2DOrigin);
return geometry2DTo1DConverter.Convert(xCoordinate);
default:
return null;
}
}
public ModelParametersForPLLines ModelParametersForPLLines
{
get { return modelParametersForPLLines; }
set { modelParametersForPLLines = value; }
}
public bool IsHydraulicShortcut
{
get { return isHydraulicShortcut; }
set { isHydraulicShortcut = value; }
}
///
/// Check if enough soil geometry data available
///
private void ThrowIfInsufficientSoilGeometryData()
{
bool hasNoGeometry1DData = (SoilGeometryType == SoilGeometryType.SoilGeometry1D) && SoilProfile == null;
bool hasNoGeometry2DData = (SoilGeometryType == SoilGeometryType.SoilGeometry2D) && (SoilGeometry2DName == null || DikeEmbankmentMaterial == null || SoilBaseDB == null);
if (hasNoGeometry1DData && hasNoGeometry2DData)
{
throw new PLLinesCreatorException("PLLinesCreator contains not enough soil geometry information (SoilProfile, SoilGeometry2DName, dikeEmbankmentMaterial or soilBase)");
}
}
///
/// Create PL2 (is pl line for penetration)
///
///
private PLLine CreatePLLine2ByExpertKnowledge(double penetrationLength, double? headInPLLine2)
{
PLLine plLine = null;
if (penetrationLength < 0.0)
{
throw new PLLinesCreatorException("Negative penetration length.");
}
if ((penetrationLength.AlmostEquals(0.0)) || (headInPLLine2 == null))
{
// No penetration, or no Head Pl2 defined, so empty pl-line will be returned
plLine = new PLLine();
}
else
{
ThrowIfInsufficientSoilGeometryData();
ThrowIfNoSurfaceLine();
ThrowIfSurfaceLineContainsNoPoints();
switch (this.SoilGeometryType)
{
case SoilGeometryType.SoilGeometry1D:
plLine = CreatePlLine2ByExpertKnowledgeFor1DGeometry(penetrationLength, headInPLLine2);
break;
case SoilGeometryType.SoilGeometry2D:
plLine = CreatePlLine2ByExpertKnowledgeFor2DGeometry(penetrationLength, headInPLLine2);
break;
}
}
return plLine;
}
///
/// Create PL2 (is pl line for penetration) for 2d-geometry
///
/// Length of the penetration.
/// The head in pl line2.
///
/// Head PL2 not defined
private PLLine CreatePlLine2ByExpertKnowledgeFor2DGeometry(double penetrationLength, double? headInPLLine2)
{
if (headInPLLine2 == null)
{
throw new PLLinesCreatorException("Head PL2 not defined");
}
PLLine plLine = new PLLine();
plLine.Points.Add(new PLLinePoint(this.surfaceLine.Geometry.Points.First().X, headInPLLine2.Value));
plLine.Points.Add(new PLLinePoint(this.surfaceLine.Geometry.Points.Last().X, headInPLLine2.Value));
return plLine;
}
///
/// Create PL2 (is pl line for penetration) for 1d-geometry
///
///
///
///
private PLLine CreatePlLine2ByExpertKnowledgeFor1DGeometry(double penetrationLength, double? headInPLLine2)
{
if (headInPLLine2 == null)
{
throw new PLLinesCreatorException("Head PL2 not defined");
}
PLLine plLine = new PLLine();
if (SoilProfile != null)
{
IList aquiferLayers = this.SoilProfile.GetAquiferLayers();
if (aquiferLayers.Count == 0)
{
throw new PLLinesCreatorException("Soil profile contains no aquifer layers.");
}
if (penetrationLength > 0 && aquiferLayers.Count > 0)
{
IList infiltrationLayers = (from SoilLayer1D layer in this.SoilProfile.Layers
where (this.SoilProfile.InBetweenAquiferLayer == null || layer.TopLevel < this.SoilProfile.InBetweenAquiferLayer.TopLevel) &&
layer.TopLevel > this.SoilProfile.BottomAquiferLayer.TopLevel
select layer).ToList();
if (infiltrationLayers.Count > 0)
{
double separationLevel = this.SoilProfile.BottomAquiferLayer.TopLevel + penetrationLength;
if (separationLevel <= infiltrationLayers.First().TopLevel)
{
plLine.Points.Add(new PLLinePoint(this.surfaceLine.Geometry.Points.First().X, headInPLLine2.Value));
plLine.Points.Add(new PLLinePoint(this.surfaceLine.Geometry.Points.Last().X, headInPLLine2.Value));
SoilLayer1D separationLayer = (from SoilLayer1D layer in infiltrationLayers
where layer.TopLevel >= separationLevel
select layer).Last();
if (separationLevel < separationLayer.TopLevel)
{
// Split layer at separation level
var extraLayer = new SoilLayer1D();
extraLayer.Assign(separationLayer);
extraLayer.Id = this.SoilProfile.GetNewUniqueLayerId();
extraLayer.TopLevel = separationLevel;
this.SoilProfile.Layers.Insert(this.SoilProfile.Layers.IndexOf(separationLayer) + 1, extraLayer);
}
}
}
this.SoilProfile.DetermineInfiltrationLayer(penetrationLength);
}
}
return plLine;
}
///
/// Check if surfaceline contains points
///
private void ThrowIfSurfaceLineContainsNoPoints()
{
if (this.surfaceLine.Geometry.Count < 1)
{
throw new PLLinesCreatorException("Surface line contains no points.");
}
}
///
/// Check if surfaceline assigned
///
private void ThrowIfNoSurfaceLine()
{
if (this.surfaceLine == null)
{
throw new PLLinesCreatorException("PLLinesCreator contains no surface line.");
}
}
///
/// Create PL3 (is phreatic level)
///
///
private PLLine CreatePLLine3ByExpertKnowledge(double headValue, double dampingFactor, double slopeGradient)
{
return CreatePLLine3Or4ByExpertKnowledge(headValue, dampingFactor, PLLineType.PL3, slopeGradient);
}
///
/// Create PL4 (is phreatic level)
///
///
private PLLine CreatePLLine4ByExpertKnowledge(double headValue, double dampingFactor, double slopeGradient)
{
return CreatePLLine3Or4ByExpertKnowledge(headValue, dampingFactor, PLLineType.PL4, slopeGradient);
}
private PLLine CreatePLLine3Or4ByExpertKnowledge(double headValue, double dampingFactor, PLLineType plLineType, double slopeGradient)
{
// Offset to solve issue MWDAM-557
const double offset = 0.01;
PLLine plLine = new PLLine();
ThrowIfInsufficientSoilGeometryData();
ThrowIfNoSurfaceLine();
ThrowIfSurfaceLineContainsNoPoints();
// Soilprofile is used to check if there is really an aquifer below toe of dike at river.
// The assumption is made that if there is an aquifer at the riverside and at the the polderside,
// that there is a connection between these aquifers.
// In the uplift calculation there will also be a check on the existence of an aquifer.
SoilProfile1D actualSoilProfile = GetSoilProfileBelowPoint(this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver).X + offset);
if (this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver) == null)
{
throw new PLLinesCreatorException("Characteristic point \"dike toe at river\" is not defined.");
}
if (dampingFactor < 0.0)
{
throw new PLLinesCreatorException("Damping factor < 0.0");
}
SoilLayer1D relevantAquiferLayer = GetRelevantAquiferLayer(plLineType, actualSoilProfile);
if (relevantAquiferLayer != null)
{
double referenceLevel = (this.HeadInPLLine2 != null) ? this.HeadInPLLine2.Value : this.waterLevelPolder;
double headAtPolderDikeToe = headValue - Math.Max(0, dampingFactor * (headValue - referenceLevel));
plLine.Points.Add(new PLLinePoint(this.surfaceLine.Geometry.Points.First().X, headValue));
plLine.Points.Add(new PLLinePoint(this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver).X, headValue));
plLine.Points.Add(new PLLinePoint(this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X, headAtPolderDikeToe));
// Now continue PLline to the end with a slope of slopegradient
AddTailOfPl3OrPl4WithSlopeGradient(slopeGradient, plLine);
if (isAdjustPL3AndPL4SoNoUpliftWillOccurEnabled)
{
AdjustLineAccordingToTRWUplift(plLine, plLineType, slopeGradient);
}
EnsureDescendingLine(plLine);
RemoveRedundantPoints(plLine);
}
return plLine;
}
///
/// Continue PLline to the end with a slope of slopegradient
/// If PLLine is lower then polderlevel, continue with polderlevel
///
/// The slope gradient.
/// The pl line.
private void AddTailOfPl3OrPl4WithSlopeGradient(double slopeGradient, PLLine plLine)
{
if (plLine.Points.Last().Z <= this.WaterLevelPolder)
{
// the PL line is already at WaterLevelPolder, so countinue it horizontally
plLine.Points.Add(new PLLinePoint(this.surfaceLine.Geometry.Points.Last().X, this.WaterLevelPolder));
}
else
{
double lengthFromLastPlPointToEnd = Math.Abs(this.surfaceLine.Geometry.Points.Last().X - plLine.Points.Last().X);
double headAtLastPlPoint = plLine.Points.Last().Z;
double headAtEnd = plLine.Points.Last().Z - slopeGradient * lengthFromLastPlPointToEnd;
Line waterLevelPolderLine =
new Line(new GeometryPoint(this.surfaceLine.Geometry.Points.First().X, 0, this.WaterLevelPolder),
new GeometryPoint(this.surfaceLine.Geometry.Points.Last().X, 0, this.WaterLevelPolder));
Line slopeLine =
new Line(
new GeometryPoint(this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X, 0, headAtLastPlPoint),
new GeometryPoint(this.surfaceLine.Geometry.Points.Last().X, 0, headAtEnd));
GeometryPoint intersectionPoint = new GeometryPoint();
if (waterLevelPolderLine.IntersectsZ(slopeLine, out intersectionPoint))
{
plLine.Points.Add(new PLLinePoint(intersectionPoint.X, this.WaterLevelPolder));
plLine.Points.Add(new PLLinePoint(this.surfaceLine.Geometry.Points.Last().X, this.WaterLevelPolder));
}
else
{
plLine.Points.Add(new PLLinePoint(this.surfaceLine.Geometry.Points.Last().X, headAtEnd));
}
}
}
private static SoilLayer1D GetRelevantAquiferLayer(PLLineType type, SoilProfile1D actualSoilProfile)
{
IList aquiferLayers = actualSoilProfile.GetAquiferLayers();
// Note : if no aquifers at all, always throw a message
if (aquiferLayers.Count == 0)
{
string message = "Soil profile (" + actualSoilProfile.Name + ") contains no aquifer layers at all.";
throw new PLLinesCreatorException(message);
}
SoilLayer1D relevantAquiferLayer = null;
switch (type)
{
case PLLineType.PL3:
relevantAquiferLayer = actualSoilProfile.BottomAquiferLayer; // Pleistocene
break;
case PLLineType.PL4:
relevantAquiferLayer = actualSoilProfile.InBetweenAquiferLayer; // Intermediate sand layer
break;
default:
throw new PLLinesCreatorException(String.Format("Invalid PL line type:{0} for creation of PL Line 3 or 4", type));
}
// Note : DAM already handles a missing (null) InBetweenAquiferLayer itself so do not throw for that.
if (relevantAquiferLayer == null && type == PLLineType.PL3)
{
string message = "Soil profile (" + actualSoilProfile.Name + ") contains no relevant aquifer layer.";
throw new PLLinesCreatorException(message);
}
return relevantAquiferLayer;
}
///
/// Remove redundant points (i.e. lying on a line between its both neighbours)
///
///
///
private static void RemoveRedundantPoints(PLLine plLine)
{
const double Tolerance = 0.001;
for (int pointIndex = plLine.Points.Count - 2; pointIndex > 0; pointIndex--)
{
PLLinePoint plPointPrev = plLine.Points[pointIndex - 1];
PLLinePoint plPoint = plLine.Points[pointIndex];
PLLinePoint plPointNext = plLine.Points[pointIndex + 1];
if (Math.Abs((plPoint.Z - plPointPrev.Z) / (plPoint.X - plPointPrev.X) - (plPointNext.Z - plPoint.Z) / (plPointNext.X - plPoint.X)) < Tolerance)
{
plLine.Points.Remove(plPoint);
}
}
}
///
/// All points should have descending z from dike to polder
///
///
private static void EnsureDescendingLine(PLLine plLine)
{
PLLinePoint plPointPrevious = null;
foreach (PLLinePoint plPoint in plLine.Points)
{
if (plPointPrevious != null && plPoint.Z > plPointPrevious.Z)
plPoint.Z = plPointPrevious.Z;
plPointPrevious = plPoint;
}
}
///
/// Clear outputvalues for PL3 or PL4
///
///
private void ClearOutputValuesForPl3_4(PLLineType plLineType)
{
switch (plLineType)
{
case PLLineType.PL3:
pl3MinUplift = Double.MaxValue;
pl3HeadAdjusted = Double.MaxValue;
pl3LocationXMinUplift = Double.MaxValue;
break;
case PLLineType.PL4:
pl4MinUplift = Double.MaxValue;
pl4HeadAdjusted = Double.MaxValue;
pl4LocationXMinUplift = Double.MaxValue;
break;
}
}
///
/// Determine outputvalues for PL3 or PL4 for minimal upliftfactor
///
/// Type of the pl line.
/// The x coordinate.
/// The uplift factor.
/// The head value.
private void UpdateOutputValuesForPl3_4(PLLineType plLineType, double xCoordinate, double upliftFactor, double headValue)
{
switch (plLineType)
{
case PLLineType.PL3:
if (upliftFactor < pl3MinUplift)
{
pl3MinUplift = upliftFactor;
pl3HeadAdjusted = headValue;
pl3LocationXMinUplift = xCoordinate;
}
break;
case PLLineType.PL4:
if (upliftFactor < pl4MinUplift)
{
pl4MinUplift = upliftFactor;
pl4HeadAdjusted = headValue;
pl4LocationXMinUplift = xCoordinate;
}
break;
}
}
///
/// Finds the index of point with X coordinate.
///
/// The pl line.
/// The x coordinate.
///
private int FindIndexOfPointInPLLineWithXCoordinate(PLLine plLine, double x)
{
for (int pointIndex = plLine.Points.Count - 1; pointIndex > 0; pointIndex--)
{
PLLinePoint plPoint = plLine.Points[pointIndex];
if (plPoint.X.AlmostEquals(x, GeometryPoint.Precision))
{
return pointIndex;
}
}
return -1;
}
///
/// Removes the index of all points of pl line after the start index.
///
/// The pl line.
/// The start index (this point will not be removed).
private void RemoveAllPointsOfPlLineAfterIndex(PLLine plLine, int startIndex)
{
for (int pointIndex = plLine.Points.Count - 1; pointIndex > startIndex; pointIndex--)
{
plLine.Points.RemoveAt(pointIndex);
}
}
///
/// Determines whether a x coordinate is in ditch (excluding first and last point of ditch).
///
/// The x coordinate.
///
private bool IsXCoordinateInDitch(double xCoordinate)
{
if (this.surfaceLine.HasDitch())
{
return (xCoordinate > this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide).X &&
xCoordinate < this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchPolderSide).X);
}
else
{
// if ditch does not exist the x coordinate is clearly not in the ditch
return false;
}
}
///
/// Determines the thickness of aquitard at characteristic point.
/// This is the distance betweeen the surfacelevel and the top of the first aquifer
///
/// Type of the characteristic point.
///
private double? DetermineThicknessAquitardAtCharacteristicPoint(CharacteristicPointType characteristicPointType)
{
if (!surfaceLine.HasAnnotation(characteristicPointType))
{
return null;
}
else
{
var characteristicGeometryPoint = surfaceLine.CharacteristicPoints.GetGeometryPoint(characteristicPointType);
SoilProfile1D actualSoilProfile = GetSoilProfileBelowPoint(characteristicGeometryPoint.X);
double bottomLevel = actualSoilProfile.GetTopLevelOfHighestAquifer();
double topLevel = characteristicGeometryPoint.Z;
return topLevel - bottomLevel;
}
}
///
/// The correction is applied as follows:
/// Check every point from DikeToeAtPolder to SurfaceLevelInside (from left to right) for uplift.
/// If uplift occurs, then correct PL3/PL4 value, so uplift does not occur.
/// All points in PL3 from this point to DikeToeAtRiver should be removed.
/// The PL3 continues from this point on with the specified slopegradient until polderlevel.
/// Make sure PL3 is always descending from left to right.
///
/// A better implementation (not implemented yet) would be:
/// Correct plline 3 or 4 for uplift according to
/// TRW (Technisch Rapport Waterspanningen bij dijken) par. b1.3.4 "Stijghoogte in het eerste watervoerende pakket"
/// - Adjust PL3/4 for all surface points from end of profile to toe of dike, so no uplift will occur in that surface point
/// - From the point, closest to the dike, (firstAdjustedPLPoint) where this correction has been made the following has to be done
/// * PL3/4 will continue horizontally from firstAdjustedPLPoint over a distance L = 2* d (d is height all layers above the aquifer)
/// * The the PL3/4 will go down in a slope of 1:50 to the PolderLevel
/// PL3/4-----
/// \___________ L = 2 * d
/// \
/// \__________
///
///
///
///
private void AdjustLineAccordingToTRWUplift(PLLine plLine, PLLineType plLineType, double slopeGradient)
{
const double cTolerancePoint = 0.0001;
ClearOutputValuesForPl3_4(plLineType);
var upliftCalculator = new UpliftCalculator
{
IsUseOvenDryUnitWeight = IsUseOvenDryUnitWeight,
UnitWeightSoilEmbankment = (this.DikeEmbankmentMaterial == null) ? (double?)null : this.DikeEmbankmentMaterial.AbovePhreaticLevel
};
GeometryPoint startSurfacePoint = this.surfaceLine.GetDikeToeInward();
int indexOfFixedPlPoint = FindIndexOfPointInPLLineWithXCoordinate(plLine, this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver).X);
if (indexOfFixedPlPoint < 0)
{
throw new PLLinesCreatorException("Could not find fixed point for PLLine");
}
// Adjust PL3/4 for all surface points from toe of dike to end of profile to, so no uplift will occur in that surface point
IEnumerable relevantSurfacePointsList = from GeometryPoint point in this.surfaceLine.Geometry.Points
where point.X >= startSurfacePoint.X
orderby point.X ascending
select point;
// Adjustment will only be applied if the value to adjust to is smaller than the previous adjusted value (to avoid that the PL3/PL4 will be adjusted to the
// polderside of the ditch i.s.o. the dikeside of the ditch.
// So we remember which was the last adjusted value in lastAdjustedHeadOfPlLine.
double lastAdjustedHeadOfPlLine = Double.MaxValue;
bool isRemoveAllPlPointsBackToDikeToeAtRiver;
bool isSkipAdjustmentInDitch;
DetermineHowToActDueToDitch(out isSkipAdjustmentInDitch, out isRemoveAllPlPointsBackToDikeToeAtRiver);
foreach (GeometryPoint surfacePoint in relevantSurfacePointsList)
{
if (IsXCoordinateInDitch(surfacePoint.X) && isSkipAdjustmentInDitch)
{
continue;
}
ConfigureUpliftCalculator(plLineType, upliftCalculator, surfacePoint);
SoilLayer1D relevantSandLayer = GetRelevantAquiferLayer(plLineType, upliftCalculator.SoilProfile);
// When plLineType is PL4 it is possible that no relevant aquifer layer is found, then the adjustment is not necessary
if (relevantSandLayer != null)
{
double aquiferTopLayer = Math.Min(surfacePoint.Z, relevantSandLayer.TopLevel);
upliftCalculator.TopOfLayerToBeEvaluated = aquiferTopLayer;
double headOfPlLine = plLine.ZFromX(surfacePoint.X);
double upliftFactor = upliftCalculator.CalculateUpliftFactor(headOfPlLine);
if (upliftFactor <= cUpliftFactorEquilibrium)
{
// Adjust headOfPLLine so there is equilibrium
bool hasNoCoverLayer = surfacePoint.Z <= aquiferTopLayer;
if (hasNoCoverLayer)
{
headOfPlLine = this.waterLevelPolder;
}
else
{
headOfPlLine = Math.Max(upliftCalculator.CalculateHeadOfPLLine(cUpliftFactorEquilibrium),
this.waterLevelPolder);
}
if (headOfPlLine < lastAdjustedHeadOfPlLine)
{
var plPoint = new PLLinePoint(surfacePoint.X, headOfPlLine);
// Remove all points of PLLine after fixed point
RemoveAllPointsOfPlLineAfterIndex(plLine, indexOfFixedPlPoint);
plLine.Points.Add(plPoint);
if (!isRemoveAllPlPointsBackToDikeToeAtRiver)
{
// To make sure that not all points of the PL-line back to the toe of the dike at riverside are to be removed
indexOfFixedPlPoint = plLine.Points.Count - 1;
}
AddTailOfPl3OrPl4WithSlopeGradient(slopeGradient, plLine);
lastAdjustedHeadOfPlLine = headOfPlLine;
}
}
UpdateOutputValuesForPl3_4(plLineType, surfacePoint.X, upliftFactor, headOfPlLine);
}
}
// Recheck on Uplift in case points are removed
if (isRemoveAllPlPointsBackToDikeToeAtRiver && this.SurfaceLine.HasDitch())
{
GeometryPoint startPoint = this.SurfaceLine.GetDikeToeInward();
GeometryPoint endPoint = this.SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchPolderSide);
IEnumerable recheckSurfacePointsList = from GeometryPoint point in this.surfaceLine.Geometry.Points
where point.X >= startPoint.X && point.X <= endPoint.X
orderby point.X ascending
select point;
foreach (GeometryPoint surfacePoint in recheckSurfacePointsList)
{
ConfigureUpliftCalculator(plLineType, upliftCalculator, surfacePoint);
SoilLayer1D relevantSandLayer = GetRelevantAquiferLayer(plLineType, upliftCalculator.SoilProfile);
// When plLineType is PL4 it is possible that no relevant aquifer layer is found, then the adjustment is not necessary
if (relevantSandLayer != null)
{
double aquiferTopLayer = Math.Min(surfacePoint.Z, relevantSandLayer.TopLevel);
upliftCalculator.TopOfLayerToBeEvaluated = aquiferTopLayer;
double headOfPlLine = plLine.ZFromX(surfacePoint.X);
double upliftFactor = upliftCalculator.CalculateUpliftFactor(headOfPlLine);
if (upliftFactor <= cUpliftFactorEquilibrium)
{
// Adjust headOfPLLine so there is equilibrium
headOfPlLine = Math.Max(upliftCalculator.CalculateHeadOfPLLine(cUpliftFactorEquilibrium), this.waterLevelPolder);
double currentHead = plLine.ZFromX(surfacePoint.X);
if (headOfPlLine < currentHead)
{
PLLinePoint plPoint = plLine.EnsurePointAtX(surfacePoint.X, cTolerancePoint);
plPoint.Z = headOfPlLine;
}
}
UpdateOutputValuesForPl3_4(plLineType, surfacePoint.X, upliftFactor, headOfPlLine);
}
}
}
}
///
/// Configures the uplift calculator.
///
/// Type of the pl line.
/// The uplift calculator.
/// The surface point.
private void ConfigureUpliftCalculator(PLLineType plLineType, UpliftCalculator upliftCalculator, GeometryPoint surfacePoint)
{
// Offset to solve issue MWDAM-764
const double offset = 0.01;
upliftCalculator.SurfaceLevel = surfacePoint.Z;
if (currentPL1Line != null)
{
upliftCalculator.PhreaticLevel = currentPL1Line.ZFromX(surfacePoint.X);
// set phreatic level to calculate upliftfactor
}
else
{
upliftCalculator.PhreaticLevel = this.WaterLevelPolder;
// if not PL1 created then assume phreatic level is same as waterlevel polder
}
SoilProfile1D actualSoilProfile = GetSoilProfileBelowPoint(surfacePoint.X + offset);
// At the end of a geometry, there are no layers to be found beyond that end. In that case
// get the layers just before the end of the geometry.
if (actualSoilProfile.LayerCount == 0)
{
actualSoilProfile = GetSoilProfileBelowPoint(surfacePoint.X - offset);
}
AdaptSoilProfileForSurfaceLevelAccuracy(actualSoilProfile, upliftCalculator.SurfaceLevel);
upliftCalculator.SoilProfile = actualSoilProfile;
}
///
/// Determines how to act due to ditch.
///
/// if set to true [is skip adjustment in ditch].
/// if set to true [is remove all pl points back to dike toe at river].
private void DetermineHowToActDueToDitch(out bool isSkipAdjustmentInDitch, out bool isRemoveAllPlPointsBackToDikeToeAtRiver)
{
if (this.surfaceLine.HasDitch())
{
double thicknessAquitardAtTopDitch =
DetermineThicknessAquitardAtCharacteristicPoint(CharacteristicPointType.DitchPolderSide).Value;
double thicknessAquitardAtBottomDitch =
DetermineThicknessAquitardAtCharacteristicPoint(CharacteristicPointType.BottomDitchPolderSide).Value;
double widthDitchAtTop = this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchPolderSide).X -
this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide).X;
double widthDitchAtTBottom =
this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.BottomDitchPolderSide).X -
this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.BottomDitchDikeSide).X;
if ((thicknessAquitardAtTopDitch > widthDitchAtTop) &&
(thicknessAquitardAtBottomDitch > widthDitchAtTBottom))
{
isRemoveAllPlPointsBackToDikeToeAtRiver = false;
isSkipAdjustmentInDitch = true;
}
else
{
isRemoveAllPlPointsBackToDikeToeAtRiver = true;
isSkipAdjustmentInDitch = false;
}
}
else
{
isRemoveAllPlPointsBackToDikeToeAtRiver = false;
isSkipAdjustmentInDitch = false;
}
}
///
/// Due to accuracy surfacelevel could be slightly below toplevel of toplevel soilprofile
/// Adjust toplevel of soilprofile to avoid this
///
///
///
private static void AdaptSoilProfileForSurfaceLevelAccuracy(SoilProfile1D actualSoilProfile, double surfaceLevel)
{
if (surfaceLevel > actualSoilProfile.TopLevel)
{
actualSoilProfile.Layers[0].TopLevel = surfaceLevel;
}
}
///
///
///
///
///
void CopyPointsInPLLine(ref PLLine plLine, GeometryPointString genericLine)
{
plLine.Points.Clear();
foreach (GeometryPoint point in genericLine.Points)
{
plLine.Points.Add(new PLLinePoint(point.X, point.Z));
}
}
///
///
///
///
///
private PLLines CreateAllPLLinesWithExpertKnowledge(Location location)
{
PLLines plLines = new PLLines();
foreach (PLLineType plLineType in Enum.GetValues(typeof(PLLineType)))
{
bool isPL1LineDefinedForLocation = (location != null) && (location.LocalXZPL1Line != null) && (location.LocalXZPL1Line.Points.Count > 1);
if ((plLineType == PLLineType.PL1) && isPL1LineDefinedForLocation)
{
PLLine plLine = plLines.Lines[plLineType];
CopyPointsInPLLine(ref plLine, location.LocalXZPL1Line);
}
else
{
plLines.Lines[plLineType] = CreatePLLineByExpertKnowledge(plLineType, location.DamType, location.SlopeDampingPiezometricHeightPolderSide);
}
// currentPL1Line is needed when calculating uplift reduction for PL3 and Pl4
if (plLineType == PLLineType.PL1)
{
// var plLine = plLines.Lines[plLineType];
// AdaptPL1ForNonWaterRetainingObject(ref plLine);
// plLines.Lines[plLineType] = plLine;
currentPL1Line = plLines.Lines[plLineType];
}
}
return plLines;
}
private IEnumerable FindAllPlLinePointsAtNWOLocation(PLLine pl1Line)
{
IEnumerable results = pl1Line.GetPointSegmentBetween(this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint1).X - cOffsetPhreaticLineBelowSurface,
this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint4).X + cOffsetPhreaticLineBelowSurface);
return results;
}
private void AdaptPL1ForNonWaterRetainingObject(ref PLLine pl1Line)
{
// check if nwo points are available as CharacteristicPoints
var nwo1 = this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint1);
var nwo2 = this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint2);
var nwo3 = this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint3);
var nwo4 = this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint4);
if ((nwo1 != null) && (nwo2 != null) && (nwo3 != null) && (nwo4 != null))
{
// Find all points in the Pl line that (might) coincide with the NWO
IEnumerable plPointsToBeMoved = FindAllPlLinePointsAtNWOLocation(pl1Line);
PLLinePoint nwo1Pl = null;
PLLinePoint nwo2Pl = null;
PLLinePoint nwo3Pl = null;
PLLinePoint nwo4Pl = null;
// For NWO point, determine whether a pl point has to be added
if (pl1Line.PositionXzOfPointRelatedToPLLine(nwo1) != PLLinePointPositionXzType.AbovePLLine)
{
nwo1Pl = new PLLinePoint(nwo1.X, nwo1.Z - cOffsetPhreaticLineBelowSurface);
}
if (pl1Line.PositionXzOfPointRelatedToPLLine(nwo2) != PLLinePointPositionXzType.AbovePLLine)
{
nwo2Pl = new PLLinePoint(nwo2.X, nwo2.Z - cOffsetPhreaticLineBelowSurface);
}
if (pl1Line.PositionXzOfPointRelatedToPLLine(nwo3) != PLLinePointPositionXzType.AbovePLLine)
{
nwo3Pl = new PLLinePoint(nwo3.X, nwo3.Z - cOffsetPhreaticLineBelowSurface);
}
if (pl1Line.PositionXzOfPointRelatedToPLLine(nwo4) != PLLinePointPositionXzType.AbovePLLine)
{
nwo4Pl = new PLLinePoint(nwo4.X, nwo4.Z - cOffsetPhreaticLineBelowSurface);
}
// Find the intersections of pl line and NWO and handle these
// Intersection between nwo point1 and nwo point2 only when nwo point1 is above pl line and nwo point2 is below plLine
PLLinePoint intersection1 = null;
if ((nwo1Pl == null) && (nwo2Pl != null))
{
var lineNWO = new Deltares.Geometry.Line { BeginPoint = new GeometryPoint(nwo1.X, 0, nwo1.Z), EndPoint = new GeometryPoint(nwo2.X, 0, nwo2.Z) };
var ips = pl1Line.IntersectionPointsXzWithLineXz(lineNWO);
if (ips.Count > 0)
{
intersection1 = new PLLinePoint(ips.First().X, ips.First().Z);
}
}
// Intersection between nwo point3 and nwo point4 only when nwo point3 is below pl line and nwo point4 is above plLine
PLLinePoint intersection2 = null;
if ((nwo3Pl != null) && (nwo4Pl == null))
{
var lineNWO = new Deltares.Geometry.Line { BeginPoint = new GeometryPoint(nwo3.X, 0, nwo3.Z), EndPoint = new GeometryPoint(nwo4.X, 0, nwo4.Z) };
var ips = pl1Line.IntersectionPointsXzWithLineXz(lineNWO);
if (ips.Count > 0)
{
intersection2 = new PLLinePoint(ips.Last().X, ips.Last().Z);
}
}
// Handle making the NWO empty
if ((NWOPhreaticAdaption != null) && (NWOPhreaticAdaption == PhreaticAdaptionType.MakeEmpty))
{
// for the polderside, the pl line is always allowed to be adapted. For the riverside, the pl line may only be adapted when the original waterlevel is runs through the nwo.
RemoveAllWaterFromNonWaterRetainingObject(nwo1, pl1Line, nwo1Pl, nwo2Pl, nwo3Pl, nwo4Pl, intersection1, intersection2, plPointsToBeMoved);
}
// Handle making the waterlevel horizontal in the NWO at the Riverside when needed (Polderside is already done when needed, see CreatePhreaticLineSegmentsInShoulderAndPolder.
if ((NWOPhreaticAdaption != null) && (NWOPhreaticAdaption == PhreaticAdaptionType.None))
{
// For the riverside, the pl line may only be adapted when the original waterlevel is runs through the nwo and is not already level.
if ((nwo1.X <= this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver).X) && ((intersection1 != null) || (intersection2 != null)))
{
double requiredWaterLevel;
// Check whether adaption of intersection points is needed
if (intersection2 == null)
{
// only intersection 1 avaialable so add intersection 2
// first see if nwo3/4 intersects, if not try nwo2/3. If still no intersection found valid level not possible, raise error
MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection1(nwo2, nwo3, nwo4, pl1Line, nwo3Pl, nwo4Pl, intersection1);
requiredWaterLevel = intersection1.Z;
}
else
{
if (intersection1 == null)
{
// only intersection 2 avaialable so add intersection 1
// first see if nwo1/2 intersects, if not try nwo2/3. If still no intersection found valid level not possible, raise error
MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection2(nwo1, nwo2, nwo3, pl1Line, nwo1Pl, nwo2Pl, nwo3Pl, intersection2);
requiredWaterLevel = intersection2.Z;
}
else
{
// intersection 1 and intersection 2 available. Only act when levels were different.
requiredWaterLevel = Math.Min(intersection1.Z, intersection2.Z);
if ((Math.Abs(intersection1.Z - intersection2.Z) > GeometryPoint.Precision))
{
if (intersection1.Z < intersection2.Z)
{
// make level in NWO intersection1.Z
MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection1And2(nwo2, nwo3, nwo4, pl1Line, nwo3Pl, intersection1, intersection2);
}
else
{
// make level in NWO intersection2.Z
MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection2And1(nwo1, nwo2, nwo3, pl1Line, nwo2Pl, intersection1, intersection2);
}
}
}
}
// Move all the points in the pl line itself that need to be moved to the horizontal proper level.
foreach (var plLinePoint in plPointsToBeMoved)
{
plLinePoint.Z = requiredWaterLevel;
}
pl1Line.DeleteCoinsidingPoints(GeometryPoint.Precision);
}
}
}
}
private void MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection2And1(GeometryPoint nwo1, GeometryPoint nwo2, GeometryPoint nwo3, PLLine pl1Line, PLLinePoint nwo2Pl, PLLinePoint intersection1, PLLinePoint intersection2)
{
var lineNWO = new Deltares.Geometry.Line { BeginPoint = new GeometryPoint(nwo1.X, 0, nwo1.Z), EndPoint = new GeometryPoint(nwo2.X, 0, nwo2.Z) };
var linePL = new Deltares.Geometry.Line { BeginPoint = new GeometryPoint(nwo1.X, 0, intersection2.Z), EndPoint = new GeometryPoint(intersection2.X, 0, intersection2.Z) };
var isp = new GeometryPoint();
if (LineHelper.GetStrictIntersectionPoint(lineNWO, linePL, ref isp))
{
var newP1 = pl1Line.EnsurePointAtX(intersection2.X, GeometryPoint.Precision);
newP1.Z = intersection2.Z;
var newP2 = pl1Line.EnsurePointAtX(isp.X - cOffsetPhreaticLineBelowSurface, GeometryPoint.Precision);
newP2.Z = intersection2.Z;
var newP3 = pl1Line.EnsurePointAtX(intersection1.X - cOffsetPhreaticLineBelowSurface, GeometryPoint.Precision);
newP3.Z = intersection1.Z;
}
else
{
var lineNWOb = new Deltares.Geometry.Line { BeginPoint = new GeometryPoint(nwo2.X, 0, nwo2.Z), EndPoint = new GeometryPoint(nwo3.X, 0, nwo3.Z) };
if (LineHelper.GetStrictIntersectionPoint(lineNWOb, linePL, ref isp))
{
var newP1 = pl1Line.EnsurePointAtX(intersection2.X, GeometryPoint.Precision);
newP1.Z = intersection2.Z;
var newP2 = pl1Line.EnsurePointAtX(isp.X - cOffsetPhreaticLineBelowSurface, GeometryPoint.Precision);
newP2.Z = intersection2.Z;
var newP3 = pl1Line.EnsurePointAtX(nwo2Pl.X, GeometryPoint.Precision);
newP3.Z = nwo2Pl.Z;
if (nwo2Pl.X > intersection1.X)
{
var newP4 = pl1Line.EnsurePointAtX(intersection1.X - cOffsetPhreaticLineBelowSurface, GeometryPoint.Precision);
newP4.Z = intersection1.Z;
}
}
else
{
throw new PLLinesCreatorException("Could not create the intersectionsection points between NWO and Phreatic line to create horizontal level.");
}
}
}
private void MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection1And2(GeometryPoint nwo2, GeometryPoint nwo3, GeometryPoint nwo4, PLLine pl1Line, PLLinePoint nwo3Pl, PLLinePoint intersection1, PLLinePoint intersection2)
{
var lineNWO = new Deltares.Geometry.Line { BeginPoint = new GeometryPoint(nwo3.X, 0, nwo3.Z), EndPoint = new GeometryPoint(nwo4.X, 0, nwo4.Z) };
var linePL = new Deltares.Geometry.Line { BeginPoint = new GeometryPoint(intersection1.X, 0, intersection1.Z), EndPoint = new GeometryPoint(nwo4.X, 0, intersection1.Z) };
var isp = new GeometryPoint();
if (LineHelper.GetStrictIntersectionPoint(lineNWO, linePL, ref isp))
{
var newP1 = pl1Line.EnsurePointAtX(intersection1.X, GeometryPoint.Precision);
newP1.Z = intersection1.Z;
var newP2 = pl1Line.EnsurePointAtX(isp.X + cOffsetPhreaticLineBelowSurface, GeometryPoint.Precision);
newP2.Z = intersection1.Z;
var newP3 = pl1Line.EnsurePointAtX(intersection2.X + cOffsetPhreaticLineBelowSurface, GeometryPoint.Precision);
newP3.Z = intersection2.Z;
}
else
{
var lineNWOb = new Deltares.Geometry.Line { BeginPoint = new GeometryPoint(nwo2.X, 0, nwo2.Z), EndPoint = new GeometryPoint(nwo3.X, 0, nwo3.Z) };
if (LineHelper.GetStrictIntersectionPoint(lineNWOb, linePL, ref isp))
{
var newP1 = pl1Line.EnsurePointAtX(intersection1.X, GeometryPoint.Precision);
newP1.Z = intersection1.Z;
var newP2 = pl1Line.EnsurePointAtX(isp.X + cOffsetPhreaticLineBelowSurface, GeometryPoint.Precision);
newP2.Z = intersection1.Z;
var newP3 = pl1Line.EnsurePointAtX(nwo3Pl.X, GeometryPoint.Precision);
newP3.Z = nwo3Pl.Z;
if (nwo3Pl.X < intersection2.X)
{
var newP4 = pl1Line.EnsurePointAtX(intersection2.X + cOffsetPhreaticLineBelowSurface, GeometryPoint.Precision);
newP4.Z = intersection2.Z;
}
}
else
{
throw new PLLinesCreatorException("Could not create the intersectionsection points between NWO and Phreatic line to create horizontal level.");
}
}
}
private void RemoveAllWaterFromNonWaterRetainingObject(GeometryPoint nwo1, PLLine pl1Line, PLLinePoint nwo1Pl, PLLinePoint nwo2Pl, PLLinePoint nwo3Pl, PLLinePoint nwo4Pl, PLLinePoint intersection1, PLLinePoint intersection2, IEnumerable plPointsToBeMoved)
{
if ((nwo1.X >= this.surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X) ||
((intersection1 != null) && (intersection2 != null)))
{
// Move all the points in the pl line itself that need to be moved to below the surfaceline.
MoveSelectedPLLinePointsBelowSurfaceLine(plPointsToBeMoved);
// now add all extra points to the pl line
if (nwo1Pl != null)
{
var newP = pl1Line.EnsurePointAtX(nwo1Pl.X, GeometryPoint.Precision);
newP.Z = nwo1Pl.Z;
}
if (nwo2Pl != null)
{
var newP = pl1Line.EnsurePointAtX(nwo2Pl.X, GeometryPoint.Precision);
newP.Z = nwo2Pl.Z;
}
if (nwo3Pl != null)
{
var newP = pl1Line.EnsurePointAtX(nwo3Pl.X, GeometryPoint.Precision);
newP.Z = nwo3Pl.Z;
}
if (nwo4Pl != null)
{
var newP = pl1Line.EnsurePointAtX(nwo4Pl.X, GeometryPoint.Precision);
newP.Z = nwo4Pl.Z;
}
// Note: for intersection points, apply offset in X direction not in Z.
if (intersection1 != null)
{
var newP = pl1Line.EnsurePointAtX(intersection1.X - cOffsetPhreaticLineBelowSurface, GeometryPoint.Precision);
newP.Z = intersection1.Z;
}
if (intersection2 != null)
{
var newP = pl1Line.EnsurePointAtX(intersection2.X + cOffsetPhreaticLineBelowSurface, GeometryPoint.Precision);
newP.Z = intersection2.Z;
}
pl1Line.DeleteCoinsidingPoints(GeometryPoint.Precision);
}
}
private void MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection2(GeometryPoint nwo1, GeometryPoint nwo2, GeometryPoint nwo3, PLLine pl1Line, PLLinePoint nwo1Pl, PLLinePoint nwo2Pl, PLLinePoint nwo3Pl, PLLinePoint intersection2)
{
var lineNWO = new Deltares.Geometry.Line { BeginPoint = new GeometryPoint(nwo1.X, 0, nwo1.Z), EndPoint = new GeometryPoint(nwo2.X, 0, nwo2.Z) };
var linePL = new Deltares.Geometry.Line { BeginPoint = new GeometryPoint(nwo1.X, 0, intersection2.Z), EndPoint = new GeometryPoint(intersection2.X, 0, intersection2.Z) };
var isp = new GeometryPoint();
if (LineHelper.GetStrictIntersectionPoint(lineNWO, linePL, ref isp))
{
var newP1 = pl1Line.EnsurePointAtX(intersection2.X, GeometryPoint.Precision);
newP1.Z = intersection2.Z;
var newP2 = pl1Line.EnsurePointAtX(isp.X - cOffsetPhreaticLineBelowSurface, GeometryPoint.Precision);
newP2.Z = intersection2.Z;
if (nwo1Pl != null)
{
var newP3 = pl1Line.EnsurePointAtX(nwo1Pl.X, GeometryPoint.Precision);
newP3.Z = nwo1Pl.Z;
}
}
else
{
var lineNWOb = new Deltares.Geometry.Line { BeginPoint = new GeometryPoint(nwo2.X, 0, nwo2.Z), EndPoint = new GeometryPoint(nwo3.X, 0, nwo3.Z) };
if (LineHelper.GetStrictIntersectionPoint(lineNWOb, linePL, ref isp))
{
var newP1 = pl1Line.EnsurePointAtX(intersection2.X, GeometryPoint.Precision);
newP1.Z = intersection2.Z;
var newP2 = pl1Line.EnsurePointAtX(isp.X - cOffsetPhreaticLineBelowSurface, GeometryPoint.Precision);
newP2.Z = intersection2.Z;
if (nwo2Pl != null)
{
var newP3 = pl1Line.EnsurePointAtX(nwo2Pl.X, GeometryPoint.Precision);
newP3.Z = nwo2Pl.Z;
if ((nwo1Pl != null) && (nwo2Pl.X > nwo1Pl.X))
{
var newP4 = pl1Line.EnsurePointAtX(nwo1Pl.X, GeometryPoint.Precision);
newP4.Z = nwo1Pl.Z;
}
}
}
else
{
throw new PLLinesCreatorException("Could not create the intersectionsection points between NWO and Phreatic line to create horizontal level.");
}
}
}
private void MakeWaterLevelHorizontalInNWOAtRiverSideUsingInterSection1(GeometryPoint nwo2, GeometryPoint nwo3, GeometryPoint nwo4, PLLine pl1Line, PLLinePoint nwo3Pl, PLLinePoint nwo4Pl, PLLinePoint intersection1)
{
var lineNWO = new Deltares.Geometry.Line { BeginPoint = new GeometryPoint(nwo3.X, 0, nwo3.Z), EndPoint = new GeometryPoint(nwo4.X, 0, nwo4.Z) };
var linePL = new Deltares.Geometry.Line { BeginPoint = new GeometryPoint(intersection1.X, 0, intersection1.Z), EndPoint = new GeometryPoint(nwo4.X, 0, intersection1.Z) };
var isp = new GeometryPoint();
if (LineHelper.GetStrictIntersectionPoint(lineNWO, linePL, ref isp))
{
var newP1 = pl1Line.EnsurePointAtX(intersection1.X, GeometryPoint.Precision);
newP1.Z = intersection1.Z;
var newP2 = pl1Line.EnsurePointAtX(isp.X + cOffsetPhreaticLineBelowSurface, GeometryPoint.Precision);
newP2.Z = intersection1.Z;
if (nwo4Pl != null)
{
var newP3 = pl1Line.EnsurePointAtX(nwo4Pl.X, GeometryPoint.Precision);
newP3.Z = nwo4Pl.Z;
}
}
else
{
var lineNWOb = new Deltares.Geometry.Line { BeginPoint = new GeometryPoint(nwo2.X, 0, nwo2.Z), EndPoint = new GeometryPoint(nwo3.X, 0, nwo3.Z) };
if (LineHelper.GetStrictIntersectionPoint(lineNWOb, linePL, ref isp))
{
var newP1 = pl1Line.EnsurePointAtX(intersection1.X, GeometryPoint.Precision);
newP1.Z = intersection1.Z;
var newP2 = pl1Line.EnsurePointAtX(isp.X + cOffsetPhreaticLineBelowSurface, GeometryPoint.Precision);
newP2.Z = intersection1.Z;
if (nwo3Pl != null)
{
var newP3 = pl1Line.EnsurePointAtX(nwo3Pl.X, GeometryPoint.Precision);
newP3.Z = nwo3Pl.Z;
if ((nwo4Pl != null) && (nwo4Pl.X > nwo3Pl.X))
{
var newP4 = pl1Line.EnsurePointAtX(nwo4Pl.X, GeometryPoint.Precision);
newP4.Z = nwo4Pl.Z;
}
}
}
else
{
throw new PLLinesCreatorException("Could not create the intersectionsection points between NWO and Phreatic line to create horizontal level.");
}
}
}
private void MoveSelectedPLLinePointsBelowSurfaceLine(IEnumerable plPointsToBeMoved)
{
foreach (var plLinePoint in plPointsToBeMoved)
{
// Determine which of these points must be moved and move them
if (this.surfaceLine.Geometry.PositionXzOfPointRelatedToExtrapolatedLine(plLinePoint) !=
RelativeXzPosition.BelowGeometricLine)
{
plLinePoint.Z = this.surfaceLine.Geometry.GetZAtX(plLinePoint.X) - cOffsetPhreaticLineBelowSurface;
}
}
}
///
///
///
///
///
private PLLines CreateAllPLLinesWithGaugesWithFallbackToExpertKnowledgeRRD(Location location)
{
var plLines = new PLLines();
foreach (PLLineType plLineType in Enum.GetValues(typeof(PLLineType)))
{
GaugePLLine gaugePLLine = (this.GaugePLLines != null) ? (this.GaugePLLines.FirstOrDefault(x => x.PLLineType == plLineType)) : null;
if (gaugePLLine != null && location != null)
{
Boolean isUseWaterLevel = ((plLineType == PLLineType.PL1) || (plLineType == PLLineType.PL2));
plLines.Lines[plLineType] = CreatePLLineFromGauges(gaugePLLine, this.Gauges, location, isUseWaterLevel);
}
else
plLines.Lines[plLineType] = CreatePLLineByExpertKnowledge(plLineType, location.DamType, location.SlopeDampingPiezometricHeightPolderSide);
// currentPL1Line is needed when calculating uplift reduction for PL3 and Pl4
if (plLineType == PLLineType.PL1)
{
currentPL1Line = plLines.Lines[plLineType];
}
}
return plLines;
}
///
/// Create all PLLines
///
///
public PLLines CreateAllPLLines(Location location)
{
PLLines plLines = new PLLines();
switch (modelParametersForPLLines.PLLineCreationMethod)
{
case PLLineCreationMethod.ExpertKnowledgeLinearInDike:
case PLLineCreationMethod.ExpertKnowledgeRRD: plLines = CreateAllPLLinesWithExpertKnowledge(location);
break;
case PLLineCreationMethod.GaugesWithFallbackToExpertKnowledgeRRD: plLines = CreateAllPLLinesWithGaugesWithFallbackToExpertKnowledgeRRD(location);
break;
}
// If PL Line2 does not exists, make it equal to PL line 4
if (!plLines.Lines[PLLineType.PL2].Exists())
{
plLines.Lines[PLLineType.PL2] = plLines.Lines[PLLineType.PL4].Clone();
}
if (!plLines.Lines[PLLineType.PL1].IsXAscending())
{
throw new PLLinesCreatorException("PLLine 1 not an X-ascending polyline");
}
return plLines;
}
///
///
///
///
///
public PLLine CreatePLLineByExpertKnowledge(PLLineType plLineType, DamType damType, double slopeGradient)
{
PLLine plLine = null;
switch (plLineType)
{
case PLLineType.PL1:
plLine = this.CreatePLLine1ByExpertKnowledge();
break;
case PLLineType.PL2:
plLine = this.CreatePLLine2ByExpertKnowledge(this.ModelParametersForPLLines.PenetrationLength, this.HeadInPLLine2);
break;
case PLLineType.PL3:
plLine = this.CreatePLLine3ByExpertKnowledge(HeadPL3TakingInAccountHydraulicShortcut, this.ModelParametersForPLLines.DampingFactorPL3, slopeGradient);
break;
case PLLineType.PL4:
plLine = this.CreatePLLine4ByExpertKnowledge(GetHeadPL4TakingInAccountHydraulicShortcut(damType), this.ModelParametersForPLLines.DampingFactorPL4, slopeGradient);
break;
}
if (plLine != null)
{
plLine.DeleteCoinsidingPoints();
}
return plLine;
}
/*
private PLLine CreatePLLineFromSensorGroup(PLLineType pLLineType, SensorLocation sensorLocation)
{
var creator = SensorPLLineCreator.CreateInstance(sensorLocation);
return creator.CreatePLLine(pLLineType);
}
*/
///
///
///
///
///
///
///
private PLLine CreatePLLineFromGauges(GaugePLLine gaugePLLine, IEnumerable gauges, Location location, Boolean useWaterLevel)
{
PLLine plLine = new PLLine();
double? gaugeWaterLevelRiver = null;
double? leftMostXAtRiverLevel = null;
int pointIndex = 0;
foreach (GaugePLLinePoint gaugePLLinePoint in gaugePLLine.Points)
{
double? localX = gaugePLLinePoint.X;
double? localZ = gaugePLLinePoint.Z;
if (gauges != null)
{
if (gaugePLLinePoint.GaugeIDX != null && gaugePLLinePoint.GaugeIDX != "")
{
Gauge gauge = (gauges.Where(x => x.Name == gaugePLLinePoint.GaugeIDX && x.Location == location)).FirstOrDefault();
if (gauge != null)
localX = gauge.LocalX;
else
{
throw new PLLinesCreatorException(String.Format("Gauge PL line {0} refers to an unknown gauge named '{1}' at X coordinate #{2}.",
gaugePLLine.PLLineType.ToString(), gaugePLLinePoint.GaugeIDX, pointIndex));
}
}
if (gaugePLLinePoint.GaugeIDZ != null && gaugePLLinePoint.GaugeIDZ != "")
{
Gauge gauge = (gauges.Where(x => x.Name == gaugePLLinePoint.GaugeIDZ && x.Location == location)).FirstOrDefault();
if (gauge != null)
{
if ((!gauge.Value.HasValue) || (gauge.Value == GaugeMissVal))
throw new PLLinesCreatorException(String.Format("Value of gauge {0} at location {1} in gauge PL line of type {2} is undefined.",
gauge.Name, location.Name, gaugePLLine.PLLineType.ToString()));
localZ = gauge.Value;
}
else
{
throw new PLLinesCreatorException(String.Format("Gauge PL line {0} refers to an unknown gauge named '{1}' at Z coordinate #{2}.",
gaugePLLine.PLLineType.ToString(), gaugePLLinePoint.GaugeIDZ, pointIndex));
}
}
}
if (!localX.HasValue)
throw new PLLinesCreatorException(String.Format("Gauge PL line {0}'s X coordinate #{1} is undefined.", gaugePLLine.PLLineType.ToString(), pointIndex));
else if (!localZ.HasValue)
throw new PLLinesCreatorException(String.Format("Gauge PL line {0}'s value #{1} is undefined.", gaugePLLine.PLLineType.ToString(), pointIndex));
else
{
if (!leftMostXAtRiverLevel.HasValue || localX > leftMostXAtRiverLevel)
{
plLine.Points.Add(new PLLinePoint(localX.Value, localZ.Value));
if (useWaterLevel)
{
// Have to account for waterlevel
if (!gaugeWaterLevelRiver.HasValue)
{
// Use first gauge as waterlevel
gaugeWaterLevelRiver = localZ.Value;
IList intersectionsAtX = this.surfaceLine.Geometry.IntersectionsXAtZ(gaugeWaterLevelRiver.Value);
if (intersectionsAtX.Count() > 0)
{
// this is where the waterlevel intersects with the dike
leftMostXAtRiverLevel = intersectionsAtX.First();
plLine.Points.Add(new PLLinePoint(leftMostXAtRiverLevel.Value, gaugeWaterLevelRiver.Value));
}
else
{
// No intersection with dike; polder is flooded
leftMostXAtRiverLevel = this.surfaceLine.Geometry.Points.Last().X;
plLine.Points.Add(new PLLinePoint(this.surfaceLine.Geometry.Points.Last().X, gaugeWaterLevelRiver.Value));
}
}
}
else
{
// In this case no intersection for this Plline with the dike will be considered
leftMostXAtRiverLevel = this.surfaceLine.Geometry.Points.First().X;
}
}
}
pointIndex++;
}
if (plLine.Points.Count() > 0)
{
plLine.Points.First().X = this.surfaceLine.Geometry.Points.First().X;
plLine.Points.Last().X = this.surfaceLine.Geometry.Points.Last().X;
}
return plLine;
}
///
/// Validate of all required characteristic points exists
///
private void ValidateRequiredCharacteristicPoints()
{
if (!SurfaceLine.HasAnnotation(CharacteristicPointType.DikeTopAtRiver))
{
throw new PLLinesCreatorException("Required characteristic point DikeTopAtRiver not found. ");
}
if (!SurfaceLine.HasAnnotation(CharacteristicPointType.DikeTopAtPolder))
{
throw new PLLinesCreatorException("Required characteristic point DikeTopAtPolder not found. ");
}
if (!SurfaceLine.HasAnnotation(CharacteristicPointType.DikeToeAtPolder))
{
throw new PLLinesCreatorException("Required characteristic point DikeToeAtPolder not found. ");
}
}
///
///
///
///
///
///
private void CreatePhreaticLineSegmentsInsideDikeForLowRiverLevel(PLLine phreaticLine, double lowWaterLevel, double highWaterLevel)
{
PLLinePoint intersectionLowWaterLevelWithDike = DetermineIntersectionBetweenWaterLevelAndDike(lowWaterLevel);
GeometryPoint pointDikeToeAtRiver = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtRiver);
switch (modelParametersForPLLines.PLLineCreationMethod)
{
case PLLineCreationMethod.ExpertKnowledgeLinearInDike:
const double cPLLineOffsetBelowSurface = 0.1;
PLLinePoint intersectionHighWaterLevelWithDike = DetermineIntersectionBetweenWaterLevelAndDike(highWaterLevel);
if (intersectionLowWaterLevelWithDike == null)
{
// Intersection is supposed to be at toe of dike when no intersection found (i.e. waterlevel below toe of dike)
intersectionLowWaterLevelWithDike = new PLLinePoint(pointDikeToeAtRiver.X, pointDikeToeAtRiver.Z);
}
// At this point both intersectionHighWaterLevelWithDike and intersectionLowWaterLevelWithDike are <> null
// or else an exception would have been thrown
PLLinePoint pointBelowHighLevel = new PLLinePoint(intersectionHighWaterLevelWithDike.X,
intersectionHighWaterLevelWithDike.Z - cPLLineOffsetBelowSurface);
//Add points below surface of dike talud riverside
foreach (GeometryPoint point in SurfaceLine.Geometry.Points.Where(
point => point.X > intersectionLowWaterLevelWithDike.X &&
point.X < intersectionHighWaterLevelWithDike.X))
{
phreaticLine.Points.Add(new PLLinePoint(point.X, point.Z - cPLLineOffsetBelowSurface));
}
// Add last point (below high riverlevel/dike section
phreaticLine.Points.Add(pointBelowHighLevel);
break;
case PLLineCreationMethod.ExpertKnowledgeRRD:
//Add points below surface of dike talud riverside until toe of dike riverside
foreach (GeometryPoint point in SurfaceLine.Geometry.Points.Where(
point => point.X > intersectionLowWaterLevelWithDike.X &&
point.X <= pointDikeToeAtRiver.X))
{
if (!surfaceLine.IsNonWaterRetainingObjectPoint(point))
{
phreaticLine.Points.Add(new PLLinePoint(point.X, point.Z - cPLLineOffsetBelowSurface));
}
}
// Add points below crest of dike
double offsetDikeTopAtRiver = this.waterLevelRiverHigh - PlLineOffsetBelowDikeTopAtRiver;
double offsetDikeTopAtPolder = this.waterLevelRiverHigh - PlLineOffsetBelowDikeTopAtPolder;
GeometryPoint pointDikeTopAtRiver = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
GeometryPoint pointDikeTopAtPolder = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
if (pointDikeTopAtRiver != null)
phreaticLine.Points.Add(new PLLinePoint(pointDikeTopAtRiver.X, offsetDikeTopAtRiver));
if (pointDikeTopAtRiver != null && pointDikeTopAtPolder != null)
{
CreateOptionalPointAtDikeCrestMiddle(phreaticLine, pointDikeTopAtRiver, pointDikeTopAtPolder);
}
if (pointDikeTopAtPolder != null)
phreaticLine.Points.Add(new PLLinePoint(pointDikeTopAtPolder.X, offsetDikeTopAtPolder));
break;
}
}
///
/// Creates the optional point at dike crest middle.
///
/// The phreatic line.
/// The point dike top at river.
/// The point dike top at polder.
private void CreateOptionalPointAtDikeCrestMiddle(PLLine phreaticLine, GeometryPoint pointDikeTopAtRiver,
GeometryPoint pointDikeTopAtPolder)
{
if (UsePlLineOffsetBelowDikeCrestMiddle.HasValue && UsePlLineOffsetBelowDikeCrestMiddle.Value && PlLineOffsetBelowDikeCrestMiddle != null)
{
var middleDikeCrestX = (pointDikeTopAtRiver.X + pointDikeTopAtPolder.X)*0.5;
var middleDikeCrestZ = waterLevelRiverHigh - PlLineOffsetBelowDikeCrestMiddle.Value;
// Check whether middleDikeCrestZ is above the surface line. If so, than use the value at the surfaceline instead.
var allZ = surfaceLine.Geometry.GetAllZatXForLine(middleDikeCrestX);
if (allZ.Count > 0)
{
var z = allZ.FirstOrDefault();
if (middleDikeCrestZ > z)
{
middleDikeCrestZ = z;
}
}
phreaticLine.Points.Add(new PLLinePoint(middleDikeCrestX, middleDikeCrestZ));
}
}
///
/// Create the phreatic line segments inside the dike
///
///
private void CreatePhreaticLineSegmentsInsideDikeForHighRiverLevel(PLLine phreaticLine)
{
double offsetDikeTopAtRiver = waterLevelRiverHigh - PlLineOffsetBelowDikeTopAtRiver;
double offsetDikeTopAtPolder = waterLevelRiverHigh - PlLineOffsetBelowDikeTopAtPolder;
// points created here are in dike
switch (modelParametersForPLLines.PLLineCreationMethod)
{
case PLLineCreationMethod.ExpertKnowledgeLinearInDike:
// No extra points below top of dike at river and top of dike at polder
break;
case PLLineCreationMethod.ExpertKnowledgeRRD:
var pointDikeTopAtRiver = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver);
if (pointDikeTopAtRiver != null)
phreaticLine.Points.Add(new PLLinePoint(pointDikeTopAtRiver.X, offsetDikeTopAtRiver));
var pointDikeTopAtPolder = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder);
if (pointDikeTopAtRiver != null && pointDikeTopAtPolder != null)
{
CreateOptionalPointAtDikeCrestMiddle(phreaticLine, pointDikeTopAtRiver, pointDikeTopAtPolder);
}
if (pointDikeTopAtPolder != null)
phreaticLine.Points.Add(new PLLinePoint(pointDikeTopAtPolder.X, offsetDikeTopAtPolder));
break;
}
}
///
///
///
///
private void CreatePhreaticLineSegmentsInShoulderAndPolder(PLLine phreaticLine)
{
PLLinePoint intersectionPolderLevelWithDike = DetermineIntersectionBetweenPolderLevelAndDike(this.WaterLevelPolder);
GeometryPoint dikeToeAtPolderPoint = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder);
double maxXCoordinateSurface = surfaceLine.Geometry.Points.Max(x => x.X);
if (modelParametersForPLLines.PLLineCreationMethod != PLLineCreationMethod.ExpertKnowledgeLinearInDike)
{
// Points created below are points starting at shoulder point to the right
GeometryPoint shoulderBaseInsidePoint = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.ShoulderBaseInside);
if (shoulderBaseInsidePoint != null)
{
double zLevel = Math.Min(phreaticLine.Points.Last().Z,
shoulderBaseInsidePoint.Z -
Math.Max(cOffsetPhreaticLineBelowSurface, PlLineOffsetBelowShoulderBaseInside));
zLevel = Math.Max(zLevel, this.WaterLevelPolder);
// Add point if it lies left of intersection of polderlevel with dike)
if ((intersectionPolderLevelWithDike == null) ||
(intersectionPolderLevelWithDike.X > shoulderBaseInsidePoint.X))
{
phreaticLine.Points.Add(new PLLinePoint(shoulderBaseInsidePoint.X, zLevel));
}
if (UsePlLineOffsetFactorBelowShoulderCrest.HasValue && UsePlLineOffsetFactorBelowShoulderCrest.Value &&
PlLineOffsetFactorBelowShoulderCrest != null && dikeToeAtPolderPoint != null)
{
GeometryPoint shoulderTopInsidePoint =
surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.ShoulderTopInside);
if (shoulderTopInsidePoint != null)
{
zLevel = dikeToeAtPolderPoint.Z + (PlLineOffsetFactorBelowShoulderCrest.Value*
(shoulderTopInsidePoint.Z - dikeToeAtPolderPoint.Z));
zLevel = Math.Min(zLevel, shoulderTopInsidePoint.Z - cOffsetPhreaticLineBelowSurface);
zLevel = Math.Min(zLevel, phreaticLine.Points.Last().Z);
zLevel = Math.Max(zLevel, this.WaterLevelPolder);
phreaticLine.Points.Add(new PLLinePoint(shoulderTopInsidePoint.X, zLevel));
}
}
}
}
GeometryPoint ditchDikeSidePoint = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide);
if (dikeToeAtPolderPoint != null)
{
double zLevel = Math.Min(phreaticLine.Points.Last().Z, dikeToeAtPolderPoint.Z - Math.Max(cOffsetPhreaticLineBelowSurface, PlLineOffsetBelowDikeToeAtPolder));
if (ditchDikeSidePoint != null)
{
if (ditchDikeSidePoint.LocationEquals(dikeToeAtPolderPoint))
{
// If DikeToeAtPolder is same as DitchDikeSide pl1 should always go to polderlevel at this point
zLevel = this.WaterLevelPolder;
}
}
zLevel = Math.Max(zLevel, this.WaterLevelPolder);
// Add point if it lies left of intersection of polderlevel with dike
if ((intersectionPolderLevelWithDike == null) || (intersectionPolderLevelWithDike.X > dikeToeAtPolderPoint.X))
{
phreaticLine.Points.Add(new PLLinePoint(dikeToeAtPolderPoint.X, zLevel));
}
}
if (intersectionPolderLevelWithDike != null)
{
phreaticLine.Points.Add(intersectionPolderLevelWithDike);
}
var isDitchPresent = (ditchDikeSidePoint != null);
var isNonWaterRetainingOjectPresent = ((NWOPhreaticAdaption != null) && surfaceLine.HasAnnotation(CharacteristicPointType.NonWaterRetainingObjectPoint1));
var adjustDitch = false;
// Handle making the waterlevel horizontal in the NWO at the Polderside when needed (For Riverside see AdaptPL1ForNonWaterRetainingObject).
var nonWaterRetainingGeometryPoint = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint1);
if (isNonWaterRetainingOjectPresent && (NWOPhreaticAdaption != PhreaticAdaptionType.Fill) &&
(nonWaterRetainingGeometryPoint.X >= surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X))
{
// if there is a ditch and it is to the left of the NWO, then only the ditch needs to be adjusted and the NWO will be correct automatically.
// if there is a ditch but it is to the right of the NWO, the NWO needs adjusting and the Ditch will be correct automatically.
// if there is no ditch then the NWO needs adjusting.
if (isDitchPresent)
{
adjustDitch = (nonWaterRetainingGeometryPoint.X >= surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchPolderSide).X);
}
if (!adjustDitch)
{
GeometryPoint nw1 = nonWaterRetainingGeometryPoint;
int surfacePointIndex = surfaceLine.Geometry.Points.IndexOf(nw1);
AdjustForDitchAndOrNonWaterRetainingObjectatPolderSide(phreaticLine, surfacePointIndex);
}
}
else
{
// No (relevant) NWO so enable handling of ditch.
// If there is a ditch but there is also a NWO to the right of it at polder side, then do not adjust the ditch. Do in all other cases.
// First see if there is a NWO.
adjustDitch = !isNonWaterRetainingOjectPresent;
if (!adjustDitch)
{
// there is a NWO, check the position
adjustDitch =
!((isDitchPresent) &&
(surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint1).X >=
surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).X) &&
(surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.NonWaterRetainingObjectPoint4).X <=
surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchDikeSide).X));
}
}
// if the NWO is not there or irrelevant and there is a ditch, then adjust it.
if ((adjustDitch) && (isDitchPresent))
{
int surfacePointIndex = surfaceLine.Geometry.Points.IndexOf(ditchDikeSidePoint);
AdjustForDitchAndOrNonWaterRetainingObjectatPolderSide(phreaticLine, surfacePointIndex);
}
else
{
// if no ditch then the PL1 will continue horizontally until the end of the profile
// another choice could be to let the PL1 go to polderlevel at the end of the profile, but that could be too optimistic for uplift.
// Discussion about this was done between Vastenburg, Knoeff and The.
// After a renewed discussion with Vastenburg, Van der Zwan and Bka, it has been decided that the PL1 should follow the
// surfaceline (with offset) until either end of surface line or polder level. Note: this is only needed when the phreatic level
// at dike toe is above polder level.
if ((!isNonWaterRetainingOjectPresent) && (!isDitchPresent))
{
if (phreaticLine.Points[phreaticLine.Points.Count - 1].Z > WaterLevelPolder)
{
AddPhreaticLineAlongSurfaceLevel(phreaticLine);
}
}
}
//Validate if endpoint surface has reached
if (phreaticLine.Points.Last().X != maxXCoordinateSurface)
{
PLLinePoint endPoint = new PLLinePoint(maxXCoordinateSurface, phreaticLine.Points.Last().Z);
phreaticLine.Points.Add(endPoint);
}
}
private void AddPhreaticLineAlongSurfaceLevel(PLLine phreaticLine)
{
// Add phreatic point at offset below every surface line point as long as depth > polder level. Else determine the
// proper position of the point at polder level (intersection) and stop.
var surfacePointIndex = surfaceLine.Geometry.Points.IndexOf(surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder)) + 1;
bool intersected = false;
for (int i = surfacePointIndex; i < surfaceLine.Geometry.Points.Count; i++)
{
double z = Math.Max(waterLevelPolder, surfaceLine.Geometry.Points[i].Z - PlLineOffsetBelowDikeToeAtPolder);
double x = surfaceLine.Geometry.Points[i].X;
if (waterLevelPolder > surfaceLine.Geometry.Points[i].Z - PlLineOffsetBelowDikeToeAtPolder)
{
// Determine intersection between would be phreatic segment and polderlevel. Add that as next point.
Line waterLevelPolderLine = new Line(new GeometryPoint(surfaceLine.Geometry.Points.First().X, 0, WaterLevelPolder),
new GeometryPoint(surfaceLine.Geometry.Points.Last().X, 0, WaterLevelPolder));
Line slopeLine = new Line(new GeometryPoint(phreaticLine.Points[phreaticLine.Points.Count - 1].X, 0, phreaticLine.Points[phreaticLine.Points.Count - 1].Z),
new GeometryPoint(surfaceLine.Geometry.Points[i].X, 0, surfaceLine.Geometry.Points[i].Z - PlLineOffsetBelowDikeToeAtPolder));
GeometryPoint intersectionPoint = new GeometryPoint();
if (waterLevelPolderLine.IntersectsZ(slopeLine, out intersectionPoint))
{
x = intersectionPoint.X;
}
intersected = true;
}
PLLinePoint point = new PLLinePoint(x, z);
phreaticLine.Points.Add(point);
if (intersected)
{
break;
}
}
}
private void AdjustForDitchAndOrNonWaterRetainingObjectatPolderSide(PLLine phreaticLine, int surfacePointIndex)
{
const double maxDouble = 99999.999;
var phreaticPolderPartialLine = new Deltares.Geometry.Line();
//#bka: hier niet langer ook starten met waterlevel als waterlevel onder bottomditch zit!
phreaticPolderPartialLine.SetBeginAndEndPoints(new GeometryPoint(phreaticLine.Points[0].X, 0, waterLevelPolder),
new GeometryPoint(maxDouble, 0, waterLevelPolder));
AddIntersectionDitchDikeSegmentPolderLevelToPhreatic(phreaticLine, surfacePointIndex, phreaticPolderPartialLine);
AddIntersectionDitchPolderSegmentPolderLevelToPhreatic(phreaticLine, phreaticPolderPartialLine);
}
///
/// Intersection between two line segments:
/// -Ditchpolder Surfaceline segment
/// -Polder waterlevel
///
///
///
private void AddIntersectionDitchPolderSegmentPolderLevelToPhreatic(PLLine phreaticLine, Deltares.Geometry.Line phreaticPolderPartialLine)
{
GeometryPoint pointDitchPolder = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DitchPolderSide);
if (pointDitchPolder != null)
{
int indexatDitchPolder = surfaceLine.Geometry.Points.IndexOf(pointDitchPolder);
var lineDitchPolderSide = new Deltares.Geometry.Line();
if (indexatDitchPolder > 1)
{
lineDitchPolderSide.SetBeginAndEndPoints(new GeometryPoint(surfaceLine.Geometry.Points[indexatDitchPolder - 1].X, 0, surfaceLine.Geometry.Points[indexatDitchPolder - 1].Z), new GeometryPoint(surfaceLine.Geometry.Points[indexatDitchPolder].X, 0, surfaceLine.Geometry.Points[indexatDitchPolder].Z));
GeometryPoint intersectDitchPolderPhreatic = new GeometryPoint();
if (LineHelper.GetStrictIntersectionPoint(lineDitchPolderSide, phreaticPolderPartialLine, ref intersectDitchPolderPhreatic))
{
phreaticLine.Points.Add(new PLLinePoint(intersectDitchPolderPhreatic.X, intersectDitchPolderPhreatic.Z));
}
else
{
phreaticLine.Points.Add(new PLLinePoint(pointDitchPolder.X, phreaticPolderPartialLine.EndPoint.Z));
}
}
else
throw new PLLinesCreatorException("DitchPolderSide should be part of a line segment");
}
}
///
/// Intersection between two line segments:
/// -DitchDike Surfaceline segment
/// -Polder waterlevel
///
///
///
///
private void AddIntersectionDitchDikeSegmentPolderLevelToPhreatic(PLLine phreaticLine, int surfacePointIndex, Deltares.Geometry.Line phreaticPolderPartialLine)
{
if (surfacePointIndex + 1 < surfaceLine.Geometry.Points.Count)
{
var lineDitchDikeSide = new Deltares.Geometry.Line();
lineDitchDikeSide.SetBeginAndEndPoints(new GeometryPoint(surfaceLine.Geometry.Points[surfacePointIndex].X, 0, surfaceLine.Geometry.Points[surfacePointIndex].Z), new GeometryPoint(surfaceLine.Geometry.Points[surfacePointIndex + 1].X, 0, surfaceLine.Geometry.Points[surfacePointIndex + 1].Z));
GeometryPoint intersectDitchDikePhreatic = new GeometryPoint();
if (LineHelper.GetStrictIntersectionPoint(lineDitchDikeSide, phreaticPolderPartialLine, ref intersectDitchDikePhreatic))
{
phreaticLine.Points.Add(new PLLinePoint(intersectDitchDikePhreatic.X, intersectDitchDikePhreatic.Z));
}
else
{
phreaticLine.Points.Add(new PLLinePoint(surfaceLine.Geometry.Points[surfacePointIndex].X, phreaticPolderPartialLine.EndPoint.Z));
}
}
else
throw new PLLinesCreatorException("Could not create DikeSementPolderLevel");
}
///
/// Intersection between two line segments:
/// - Horizontal waterlevel
/// - Surface line (until DikeTopAtRiver)
///
///
///
///
public PLLinePoint DetermineIntersectionBetweenWaterLevelAndDike(double waterLevelRiver)
{
GeometryPoint PointAtLevel = this.surfaceLine.DetermineIntersectionWithLevel(waterLevelRiver);
if (PointAtLevel != null)
{
return new PLLinePoint(PointAtLevel.X, PointAtLevel.Z);
}
else
{
return null;
}
}
private PLLinePoint DetermineIntersectionBetweenPolderLevelAndDike(double polderLevel)
{
var polderlevelLine = new Deltares.Geometry.Line();
double startXCoordinate = this.surfaceLine.Geometry.Points.OrderBy(p => p.X).First().X;
GeometryPoint pointEndOfprofile = SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.SurfaceLevelInside);
polderlevelLine.SetBeginAndEndPoints(new GeometryPoint(startXCoordinate, 0, polderLevel), new GeometryPoint(pointEndOfprofile.X, 0, polderLevel));
ThrowWhenWaterLevelAboveDike(polderLevel, SurfaceLine);
int startPosition = 0;
int endPosition = 0;
startPosition = SurfaceLine.Geometry.Points.IndexOf(SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver));
endPosition = SurfaceLine.Geometry.Points.IndexOf(SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder));
for (int surfacePointIndex = startPosition; surfacePointIndex < endPosition; surfacePointIndex++)
{
var surfaceLineSegment = new Deltares.Geometry.Line();
surfaceLineSegment.SetBeginAndEndPoints(new GeometryPoint(SurfaceLine.Geometry.Points[surfacePointIndex].X, 0, SurfaceLine.Geometry.Points[surfacePointIndex].Z), new GeometryPoint(SurfaceLine.Geometry.Points[surfacePointIndex + 1].X, 0, SurfaceLine.Geometry.Points[surfacePointIndex + 1].Z));
GeometryPoint intersectPoint = new GeometryPoint();
if (LineHelper.GetStrictIntersectionPoint(surfaceLineSegment, polderlevelLine, ref intersectPoint))
{
return new PLLinePoint(intersectPoint.X, intersectPoint.Z);
}
}
return null;
}
///
/// Detects whether the phreatic line is above the surface line between DikeTopAtRiver and DikeToeAtPolder
/// If the phreatic line is above the surface, then the phreatic line should follow the surface
///
///
///
private void ValidatePhreaticBelowDike(PLLine phreaticLine)
{
int startIndex = SurfaceLine.Geometry.Points.IndexOf(SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver));
int stopIndex = SurfaceLine.Geometry.Points.IndexOf(SurfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder));
bool StopIteration = false;
// This code can go into an endless loop
// Enter failsave to raise exception when this occurs
int cMaxIterations = Math.Max(100, phreaticLine.Points.Count * SurfaceLine.Geometry.Points.Count);
int iterationCount = 0;
while (StopIteration == false)
{
iterationCount++;
bool foundIntersect = false;
for (int phreaticPointIndex = 0; phreaticPointIndex < phreaticLine.Points.Count - 1; phreaticPointIndex++)
{
var phreaticLineSegment = new Deltares.Geometry.Line();
phreaticLineSegment.SetBeginAndEndPoints(phreaticLine.Points[phreaticPointIndex], phreaticLine.Points[phreaticPointIndex + 1]);
for (int surfacePointIndex = startIndex; surfacePointIndex < stopIndex; surfacePointIndex++)
{
var surfaceLineSegment = new Deltares.Geometry.Line();
surfaceLineSegment.SetBeginAndEndPoints(SurfaceLine.Geometry.Points[surfacePointIndex], SurfaceLine.Geometry.Points[surfacePointIndex + 1]);
GeometryPoint intersectPoint = new GeometryPoint();
if (LineHelper.GetStrictIntersectionPoint(phreaticLineSegment, surfaceLineSegment, ref intersectPoint))
{
// Prevent any adding when intersectPoint is already on Pl
if (!intersectPoint.LocationEquals(phreaticLineSegment.BeginPoint) &&
!intersectPoint.LocationEquals(phreaticLineSegment.EndPoint))
{
// Only add intersectPoint if it is higher than waterlevelPolder, because the intersection could be caused
// by the fact that the waterlevel is higher than DikeToeAtPolder
if (intersectPoint.Z > waterLevelPolder + cOffsetPhreaticLineBelowSurface)
{
var newPLLinePoint = new PLLinePoint(intersectPoint.X,
intersectPoint.Z - cOffsetPhreaticLineBelowSurface);
// Only add new point if it is more to the right than the first point of the phreatic line segment when
if (newPLLinePoint.X > phreaticLineSegment.BeginPoint.X)
{
phreaticLine.Points.Insert(phreaticPointIndex + 1, newPLLinePoint);
if (SurfaceLine.Geometry.Points[surfacePointIndex + 1].X.IsNearEqual(intersectPoint.X))
{
// phreatic point and surfaceline point are postioned above each other, so replace the phreatic point with the intersect point
phreaticLine.Points[phreaticPointIndex + 2] =
new PLLinePoint(SurfaceLine.Geometry.Points[surfacePointIndex + 1].X,
SurfaceLine.Geometry.Points[surfacePointIndex + 1].Z -
cOffsetPhreaticLineBelowSurface);
}
else
{
var newNextPLLinePoint =
new PLLinePoint(SurfaceLine.Geometry.Points[surfacePointIndex + 1].X,
SurfaceLine.Geometry.Points[surfacePointIndex + 1].Z -
cOffsetPhreaticLineBelowSurface);
if (newNextPLLinePoint.X <= phreaticLine.Points[phreaticPointIndex + 2].X)
{
if (newNextPLLinePoint.X.IsNearEqual(phreaticLine.Points[phreaticPointIndex + 2].X))
{
// If phreatic point already exist on this x-coordinate just replace it with the new point
phreaticLine.Points[phreaticPointIndex + 2] = newNextPLLinePoint;
}
else
{
phreaticLine.Points.Insert(phreaticPointIndex + 2,
newNextPLLinePoint);
}
}
}
foundIntersect = true;
}
break;
}
}
}
}
}
if (foundIntersect == false)
StopIteration = true;
if (iterationCount > cMaxIterations)
throw new PLLinesCreatorException("PLLinesCreator.ValidatePhreaticBelowDike() cannot succeed");
}
}
///
/// Create PL1 (is phreatic level)
///
///
///
private PLLine CreatePLLine1ByExpertKnowledge()
{
ValidateSurfaceLine();
ValidateRequiredCharacteristicPoints();
ThrowWhenWaterLevelPolderAboveDikeTopAtPolder();
//Create Phreatic line and add polderwater level
double startXCoordinate = this.surfaceLine.Geometry.Points.OrderBy(p => p.X).First().X;
PLLine phreaticLine = new PLLine();
phreaticLine.IsPhreatic = true;
double? waterLevel = WaterLevelToUse();
// Waterlevel is supposed to be at level of SurfaceLevelOutside when waterlevel is below SurfaceLevelOutside (is the bottom of the river)
waterLevel = this.SurfaceLine.EnsureWaterLevelIsAboveRiverBottom(waterLevel);
phreaticLine.Points.Add(new PLLinePoint(startXCoordinate, waterLevel.Value));
PLLinePoint intersectionPhreaticDike = null;
//Determine intersection Phreaticlevel - Berm (or dike)
intersectionPhreaticDike = DetermineIntersectionBetweenWaterLevelAndDike(waterLevel.Value);
if (intersectionPhreaticDike == null)
{
throw new PLLinesCreatorException("DetermineIntersectionBetweenWaterLevelAndDike failed");
}
//Add intersectioncoordinate to phreatic line list if not already in list.
if (!phreaticLine.Points.Contains(intersectionPhreaticDike))
{
phreaticLine.Points.Add(intersectionPhreaticDike);
}
//Complete phreatic line
if (this.IsUseLowWaterLevel)
{
// intersectionPhreaticDike.Z is either low riverlevel or ToeAtDike.Z, whichever is higher
// Or in human language: low riverlevel should always be a least at surfacelevel
CreatePhreaticLineSegmentsInsideDikeForLowRiverLevel(phreaticLine, Math.Max(intersectionPhreaticDike.Z, this.waterLevelRiverLow.Value), this.waterLevelRiverHigh);
}
else
{
CreatePhreaticLineSegmentsInsideDikeForHighRiverLevel(phreaticLine);
}
CreatePhreaticLineSegmentsInShoulderAndPolder(phreaticLine);
//Check if phreatic line is above
ValidatePhreaticAboveWaterLevelPolder(phreaticLine);
//Check if phreatic line is below the surface line
ValidatePhreaticBelowDike(phreaticLine);
// currentPL1Line is needed when calculating uplift reduction for PL3 and Pl4
AdaptPL1ForNonWaterRetainingObject(ref phreaticLine);
currentPL1Line = phreaticLine;
return phreaticLine;
}
///
/// Check whether a correct surface line has been imported
///
private void ValidateSurfaceLine()
{
if (SurfaceLine == null)
{
throw new PLLinesCreatorException("Surfaceline should be initialized");
}
if (SurfaceLine.Geometry.Points.Count < 2)
{
throw new PLLinesCreatorException("More surfaceline points are required");
}
}
///
/// Check if phreatic line is above waterlevel
/// If it is below, correct it
///
///
private void ValidatePhreaticAboveWaterLevelPolder(PLLine phreaticLine)
{
foreach (PLLinePoint phreaticPoint in phreaticLine.Points)
{
if (phreaticPoint.Z < waterLevelPolder)
{
phreaticPoint.Z = waterLevelPolder;
}
}
}
///
///
///
private void ThrowWhenWaterLevelAboveDike(double waterLevelRiver, SurfaceLine2 surfaceLine)
{
var dikeTopAtRiverHeight = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtRiver).Z;
if (waterLevelRiver > dikeTopAtRiverHeight)
{
throw new PLLinesCreatorException(String.Format("Phreatic water level should have an intersection with the dike at least once (phreatic line higher ({0:F2} m) than surface line ({1:F2}))", waterLevelRiver, dikeTopAtRiverHeight));
}
}
///
///
///
private void ThrowWhenWaterLevelPolderAboveDikeTopAtPolder()
{
if (waterLevelPolder > surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder).Z)
{
throw new PLLinesCreatorException(String.Format("Waterlevel ({0:0.00}) in polder higher than dike top at polder ({1:0.00}) in surfaceline {2}", waterLevelPolder, SurfaceLine.GetDikeToeInward().Z, SurfaceLine.Name));
}
}
}
}