// Copyright (C) Stichting Deltares 2018. 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.IO; using System.Linq; using System.Text; using Deltares.DamEngine.Calculators.KernelWrappers.Common; using Deltares.DamEngine.Calculators.KernelWrappers.Interfaces; using Deltares.DamEngine.Calculators.Properties; using Deltares.DamEngine.Data.Design; using Deltares.DamEngine.Data.General; using Deltares.DamEngine.Data.General.Results; using Deltares.DamEngine.Data.Geometry; using Deltares.DamEngine.Data.Geotechnics; using Deltares.DamEngine.Data.Standard.Calculation; using Deltares.DamEngine.Data.Standard.Logging; namespace Deltares.DamEngine.Calculators.DikesDesign { /// /// The Dam Engine design calculator /// public class DesignCalculator { private const double CMinimumShoulderElevation = 0.5; private const double CMinimumShoulderExtraElevation = 0.05; private const double CMinimumShoulderWidth = 2; private const double CMinimumShoulderExtraWidth = 1; private const double CToleranceShoulderChanges = 0.001; /// /// Performs the design calculation /// /// The dam project data. /// public List Execute(DamProjectData damProjectData) { IKernelWrapper kernelWrapper = KernelWrapperHelper.CreateKernelWrapper(damProjectData.DamProjectCalculationSpecification.CurrentSpecification); if (kernelWrapper == null) { throw new NotImplementedException(Resources.DesignCalculatorKernelNotImplemented); } damProjectData.DesignCalculations = new List(); var calculationMessages = new List(); for (int locationIndex = 0; locationIndex < damProjectData.Dike.Locations.Count; locationIndex++) { var location = damProjectData.Dike.Locations[locationIndex]; for (int subSoilScenarioIndex = 0; subSoilScenarioIndex < location.Segment.SoilProfileProbabilities.Count; subSoilScenarioIndex++) { var soiProfileProbability = location.Segment.SoilProfileProbabilities[subSoilScenarioIndex]; for (int designScenarioIndex = 0; designScenarioIndex < location.Scenarios.Count; designScenarioIndex++) { // Prepare input DesignScenario designScenario = location.Scenarios[designScenarioIndex]; var damKernelInput = new DamKernelInput(); var projectPath = damProjectData.ProjectPath != "" ? damProjectData.ProjectPath : Directory.GetCurrentDirectory(); damKernelInput.ProjectDir = projectPath; damKernelInput.CalculationDir = Path.Combine(projectPath, damProjectData.CalculationMap); damKernelInput.Location = location; damKernelInput.SubSoilScenario = soiProfileProbability; //damKernelInput.DesignScenario = location.Scenarios[designScenarioIndex]; damKernelInput.DamFailureMechanismeCalculationSpecification = damProjectData.DamProjectCalculationSpecification.CurrentSpecification; damKernelInput.RiverLevelHigh = designScenario.RiverLevel; damKernelInput.RiverLevelLow = designScenario.RiverLevelLow; damKernelInput.FilenamePrefix = String.Format("Loc({0})_Sce({1})", designScenario.Location.Name, designScenario.LocationScenarioID); AnalysisType analysisType = DamProjectCalculationSpecification.SelectedAnalysisType; SynchronizeScenarioDataWithLocationData(designScenario, location); IKernelDataInput kernelDataInput; IKernelDataOutput kernelDataOutput; PrepareResult prepareResult = kernelWrapper.Prepare(damKernelInput, 0, out kernelDataInput, out kernelDataOutput); // Sometimes the kernelDataInput is not created (p.e when soilprofileprobablility is meant for // stability where Piping calc is wanted). In that case, do nothing but just skip. if (prepareResult == PrepareResult.Successful) { switch (analysisType) { case AnalysisType.AdaptGeometry: PerformDesignCalculationShoulderIterativePerPoint(kernelWrapper, kernelDataInput, kernelDataOutput, damKernelInput, designScenario, calculationMessages, damProjectData.DesignCalculations); break; case AnalysisType.NoAdaption: PerformSingleCalculation(kernelWrapper, kernelDataInput, kernelDataOutput, damKernelInput, designScenario, calculationMessages, damProjectData.DesignCalculations); break; } } else { if (prepareResult == PrepareResult.NotRelevant) { calculationMessages.Add(new LogMessage(LogMessageType.Info, null, string.Format(Resources.DesignCalculatorIrrelevant, location.Name, soiProfileProbability.ToString()))); } if (prepareResult == PrepareResult.Failed) { calculationMessages.Add(new LogMessage(LogMessageType.Error, null, string.Format(Resources.DesignCalculatorPrepareError, location.Name, soiProfileProbability.ToString()))); } } } } } return calculationMessages; } /// /// Synchronizes the scenario data with location data. /// Note that scenario data is leading when available. /// /// The scenario. /// The location. private void SynchronizeScenarioDataWithLocationData(DesignScenario designScenario, Location location) { if (designScenario.PlLineOffsetBelowDikeToeAtPolder.HasValue) { location.PlLineOffsetBelowDikeToeAtPolder = designScenario.PlLineOffsetBelowDikeToeAtPolder.Value; } if (designScenario.PlLineOffsetBelowDikeTopAtPolder.HasValue) { location.PlLineOffsetBelowDikeTopAtPolder = designScenario.PlLineOffsetBelowDikeTopAtPolder.Value; } if (designScenario.PlLineOffsetBelowDikeTopAtRiver.HasValue) { location.PlLineOffsetBelowDikeTopAtRiver = designScenario.PlLineOffsetBelowDikeTopAtRiver.Value; } if (designScenario.PlLineOffsetBelowShoulderBaseInside.HasValue) { location.PlLineOffsetBelowShoulderBaseInside = designScenario.PlLineOffsetBelowShoulderBaseInside.Value; } if (designScenario.PlLineOffsetBelowDikeCrestMiddle.HasValue) { location.PlLineOffsetBelowDikeCrestMiddle = designScenario.PlLineOffsetBelowDikeCrestMiddle; } if (designScenario.PlLineOffsetFactorBelowShoulderCrest.HasValue) { location.PlLineOffsetFactorBelowShoulderCrest = designScenario.PlLineOffsetFactorBelowShoulderCrest; } if (designScenario.UsePlLineOffsetBelowDikeCrestMiddle.HasValue) { location.UsePlLineOffsetBelowDikeCrestMiddle = designScenario.UsePlLineOffsetBelowDikeCrestMiddle; } if (designScenario.UsePlLineOffsetFactorBelowShoulderCrest.HasValue) { location.UsePlLineOffsetFactorBelowShoulderCrest = designScenario.UsePlLineOffsetFactorBelowShoulderCrest; } if (designScenario.HeadPl3.HasValue) { location.HeadPl3 = designScenario.HeadPl3.Value; } if (designScenario.HeadPl4.HasValue) { location.HeadPl4 = designScenario.HeadPl4.Value; } // Synchronize piping design parameters if (designScenario.UpliftCriterionPiping.HasValue) { location.UpliftCriterionPiping = designScenario.UpliftCriterionPiping.Value; } if (designScenario.RequiredSafetyFactorPiping.HasValue) { location.ModelFactors.RequiredSafetyFactorPiping = designScenario.RequiredSafetyFactorPiping.Value; } // Synchronize stability design parameters if (designScenario.UpliftCriterionStability.HasValue) { location.UpliftCriterionStability = designScenario.UpliftCriterionStability.Value; } if (designScenario.RequiredSafetyFactorStabilityInnerSlope.HasValue) { location.ModelFactors.RequiredSafetyFactorStabilityInnerSlope = designScenario.RequiredSafetyFactorStabilityInnerSlope.Value; } if (designScenario.RequiredSafetyFactorStabilityOuterSlope.HasValue) { location.ModelFactors.RequiredSafetyFactorStabilityOuterSlope = designScenario.RequiredSafetyFactorStabilityOuterSlope.Value; } } private static void PerformSingleCalculation(IKernelWrapper kernelWrapper, IKernelDataInput kernelDataInput, IKernelDataOutput kernelDataOutput, DamKernelInput damKernelInput, DesignScenario designScenario, List calculationMessages, List designCalculations) { // Perform validation var designResults = new List(); List locationCalculationMessages = new List(); List validationMessages; try { int errorCount = kernelWrapper.Validate(kernelDataInput, kernelDataOutput, out validationMessages); if (errorCount > 0) { locationCalculationMessages.Add(new LogMessage(LogMessageType.Error, null, string.Format(Resources.DesignCalculatorValidationFailed, damKernelInput.Location.Name, damKernelInput.SubSoilScenario.ToString()))); locationCalculationMessages.AddRange(validationMessages); } else { // Perform calculation kernelWrapper.Execute(kernelDataInput, kernelDataOutput, out locationCalculationMessages); } // Process output calculationMessages.AddRange(locationCalculationMessages); StringBuilder sb = new StringBuilder(); foreach (var message in locationCalculationMessages) { sb.Append(message.Message + Environment.NewLine); } string resultMessage = sb.ToString(); kernelWrapper.PostProcess(damKernelInput, kernelDataOutput, designScenario, resultMessage, out designResults); } catch (Exception exception) { string resultMessage = exception.Message; kernelWrapper.PostProcess(damKernelInput, kernelDataOutput, designScenario, resultMessage, out designResults); } finally { foreach (var designResult in designResults) { designCalculations.Add(designResult); } } } /// /// Ensures that the points on the surface line are never more than cDiff (0.5) apart. /// /// /// private static IEnumerable GetCheckedSurfaceLine(IEnumerable originalLine) { const double cDiff = 0.5; var newLine = new List(); double X = originalLine.First().X; foreach (var point in originalLine) { while (point.X > X + cDiff) { var newPoint = new GeometryPoint(point) { X = X + cDiff }; if (newPoint.X > newLine.Last().X) { newPoint.Z = newLine.Last().Z + ((newPoint.X - newLine.Last().X) / (point.X - newLine.Last().X)) * (point.Z - newLine.Last().Z); newLine.Add(newPoint); } X = newPoint.X; } newLine.Add(point); } return newLine; } /// /// Calculates the maximum level for the shoulder. /// /// The surface line. /// The fraction of dike height to determine maximimum shoulder height. /// private static double CalculateMaximumShoulderLevel(SurfaceLine2 surfaceLine, double maxFractionOfDikeHeightForShoulderHeight) { var top = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeTopAtPolder).Z; var bottom = surfaceLine.CharacteristicPoints.GetGeometryPoint(CharacteristicPointType.DikeToeAtPolder).Z; if (top - bottom <= 0) { throw new DesignCalculatorException(Resources.SurfaceLineShoulderAdapterMaxShoulderHeightError); } double maxHeight = Math.Abs((top - bottom) * maxFractionOfDikeHeightForShoulderHeight); return bottom + maxHeight; } /// /// Performs the design calculation shoulder iterative per point. /// This is a design strategy used for Piping /// /// The kernel wrapper. /// The kernel data input. /// The kernel data output. /// The dam kernel input. /// The calculation messages. /// The design calculations. /// private static void PerformDesignCalculationShoulderIterativePerPoint(IKernelWrapper kernelWrapper, IKernelDataInput kernelDataInput, IKernelDataOutput kernelDataOutput, DamKernelInput damKernelInput, DesignScenario designScenario, List calculationMessages, List designCalculations) { var designResults = new List(); var location = damKernelInput.Location; var surfaceLine = location.SurfaceLine; var orgSurfaceLine = surfaceLine; try { List locationCalculationMessages = new List(); double orgShoulderLength = surfaceLine.DetermineShoulderWidth(); double orgShoulderHeight = surfaceLine.DetermineShoulderHeight(); GeometryPoint startSurfacePoint = surfaceLine.GetDikeToeInward(); IEnumerable relevantSurfacePointsList = from GeometryPoint point in surfaceLine.Geometry.Points where point.X >= startSurfacePoint.X orderby point.X ascending select point; relevantSurfacePointsList = GetCheckedSurfaceLine(relevantSurfacePointsList); double oldDesiredShoulderLength = 0.0; double desiredShoulderLength = 0.0; double desiredShoulderHeight = 0.0; double oldDesiredShoulderHeight = 0.0; int pointCount = 0; foreach (var point in relevantSurfacePointsList) { pointCount++; // Calculate the piping design at the given point. This returns the required adaption (berm length and height) if any. ShoulderDesign shoulderDesign = kernelWrapper.CalculateDesignAtPoint(damKernelInput, kernelDataInput, kernelDataOutput, point, out locationCalculationMessages); if (shoulderDesign != null) { // Piping is an issue so adapt the surfaceline for it desiredShoulderLength = shoulderDesign.ShoulderLengthFromToe; desiredShoulderLength = Math.Max(desiredShoulderLength, oldDesiredShoulderLength); oldDesiredShoulderLength = desiredShoulderLength; // shoulder height is height above surfacelevel!! desiredShoulderHeight = shoulderDesign.ShoulderHeightFromToe; desiredShoulderHeight = Math.Max(desiredShoulderHeight, oldDesiredShoulderHeight); oldDesiredShoulderHeight = desiredShoulderHeight; } } if (desiredShoulderLength > 0) { desiredShoulderLength = Math.Max(desiredShoulderLength, CMinimumShoulderWidth); } if (desiredShoulderLength > 0) { desiredShoulderHeight = Math.Max(desiredShoulderHeight, CMinimumShoulderElevation); } bool isNewShoulderSameAsOriginal = ((Math.Abs(desiredShoulderLength - orgShoulderLength) < CToleranceShoulderChanges) && (Math.Abs(desiredShoulderHeight - orgShoulderHeight) < CToleranceShoulderChanges)); SurfaceLine2 newSurfaceLine = null; if (isNewShoulderSameAsOriginal) { newSurfaceLine = surfaceLine; } else { // Adapt the surfaceline for the finally required shoulder dimensions. double maxShoulderLevel = CalculateMaximumShoulderLevel(surfaceLine, 1.0); // no limit to height of shoulder var surfaceLineShoulderAdapter = new SurfaceLineShoulderAdapter(surfaceLine, location); surfaceLineShoulderAdapter.MaxShoulderLevel = maxShoulderLevel; newSurfaceLine = surfaceLineShoulderAdapter.ConstructNewSurfaceLine(desiredShoulderLength, desiredShoulderHeight, true); } damKernelInput.Location.SurfaceLine = newSurfaceLine; kernelWrapper.Prepare(damKernelInput, 0, out kernelDataInput, out kernelDataOutput); kernelWrapper.Execute(kernelDataInput, kernelDataOutput, out locationCalculationMessages); // Process output calculationMessages.AddRange(locationCalculationMessages); StringBuilder sb = new StringBuilder(); foreach (var message in locationCalculationMessages) { sb.Append(message.Message + Environment.NewLine); } string resultMessage = sb.ToString(); kernelWrapper.PostProcess(damKernelInput, kernelDataOutput, designScenario, resultMessage, out designResults); } catch (Exception exception) { string resultMessage = exception.Message; kernelWrapper.PostProcess(damKernelInput, kernelDataOutput, designScenario, resultMessage, out designResults); } finally { damKernelInput.Location.SurfaceLine = orgSurfaceLine; foreach (var designResult in designResults) { designCalculations.Add(designResult); } } string evaluationMessage; bool designSuccessful = kernelWrapper.EvaluateDesign(damKernelInput, kernelDataInput, kernelDataOutput, out evaluationMessage); if (!designSuccessful) { throw new DesignCalculatorException(String.Format(Resources.DesignUnsuccessful, damKernelInput.Location.Name, damKernelInput.SubSoilScenario) + ", " + evaluationMessage); } } } }