using System;
using System.Collections;
using System.Collections.Generic;
using Deltares.Standard;
using Deltares.Standard.Language;
using Deltares.Standard.Logging;
// =======================================================================================================================
// Class name: TUpliftPlane
//
// Description:
//
// Copyright (c) 1997-2010 Deltares
//
// Date ID Modification
// 2008-07-01 Best Created
// =======================================================================================================================
namespace Deltares.Stability.Calculation.Inner
{
public class TUpliftPlane : TSlipPlane
{
// destructor Destroy; override;
// properties
private double FActiveDMoment = 0;
private double FActiveForce = 0;
private double FActiveForce0 = 0;
private double FActiveRMoment = 0;
private double FAngle = 0;
private double FDrivingMoment = 0;
private string FErrorNote = String.Empty;
private double FHorForce = 0;
private double FHorForce0 = 0;
private double FLoadMomentLeft = 0;
private double FLoadMomentRight = 0;
private bool FPC_Ring = false;
private double FPassiveDMoment = 0;
private double FPassiveForce = 0;
private double FPassiveForce0 = 0;
private double FPassiveRMoment = 0;
private int FPlotNumber = 0;
private double FRadLeft = 0;
private double FRadRight = 0;
private double FRemoldFactor = 0;
private double FResistingMom = 0;
private double FResistingMom0 = 0;
private bool FReverseGeom = false;
private double FSchemaFactor = 0;
private double FSoilMomentLeft = 0;
private double FSoilMomentRight = 0;
private double FTextileMoment = 0;
private double FTotDrivingLeft = 0;
private double FTotDrivingRight = 0;
private int FWaterlevelNumber = 0;
private double FWatermomentLeft = 0;
private double FWatermomentRight = 0;
private double FXActiveCentre = 0;
private double FXMidLeft = 0;
private double FXMidRight = 0;
private double FXPassiveCentre = 0;
private double FZActiveCentre = 0;
private double FZMidLeft = 0;
private double FZMidRight = 0;
private double FZPassiveCentre = 0;
private double FZTangent = 0;
public TUpliftPlane()
: base()
{
// TODO: Add any constructor code here
}
public double XActiveCentre
{
get
{
return FXActiveCentre;
}
set
{
FXActiveCentre = value;
}
}
public double ZActiveCentre
{
get
{
return FZActiveCentre;
}
set
{
FZActiveCentre = value;
}
}
public double XPassiveCentre
{
get
{
return FXPassiveCentre;
}
set
{
FXPassiveCentre = value;
}
}
public double ZPassiveCentre
{
get
{
return FZPassiveCentre;
}
set
{
FZPassiveCentre = value;
}
}
public double ActiveRadius
{
get
{
return GetActiveRadius();
}
}
public double PassiveRadius
{
get
{
return GetPassiveRadius();
}
}
public double ZTangent
{
get
{
return FZTangent;
}
set
{
FZTangent = value;
}
}
public double SafetyFactor
{
get
{
return GetStabilityFactor();
}
}
public double ReliabilityIndex
{
get
{
return GetReliabilityIndex();
}
}
public double TextileMoment
{
get
{
return FTextileMoment;
}
set
{
FTextileMoment = value;
}
}
public double ActiveDMoment
{
get
{
return FActiveDMoment;
}
set
{
FActiveDMoment = value;
}
}
public double PassiveDMoment
{
get
{
return FPassiveDMoment;
}
set
{
FPassiveDMoment = value;
}
}
public double ActiveRMoment
{
get
{
return FActiveRMoment;
}
set
{
FActiveRMoment = value;
}
}
public double PassiveRMoment
{
get
{
return FPassiveRMoment;
}
set
{
FPassiveRMoment = value;
}
}
public double ResistingMom0
{
get
{
return FResistingMom0;
}
set
{
FResistingMom0 = value;
}
}
public double ResistingMom
{
get
{
return FResistingMom;
}
set
{
FResistingMom = value;
}
}
public double DrivingMoment
{
get
{
return FDrivingMoment;
}
set
{
FDrivingMoment = value;
}
}
public double ActiveForce
{
get
{
return FActiveForce;
}
}
public double PassiveForce
{
get
{
return FPassiveForce;
}
}
public double HorForce
{
get
{
return FHorForce;
}
}
public double ActiveForce0
{
get
{
return FActiveForce0;
}
}
public double PassiveForce0
{
get
{
return FPassiveForce0;
}
}
public double HorForce0
{
get
{
return FHorForce0;
}
}
public double WatermomentLeft
{
get
{
return FWatermomentLeft;
}
set
{
FWatermomentLeft = value;
}
}
public double WatermomentRight
{
get
{
return FWatermomentRight;
}
set
{
FWatermomentRight = value;
}
}
public double LoadMomentRight
{
get
{
return FLoadMomentRight;
}
set
{
FLoadMomentRight = value;
}
}
public double LoadMomentLeft
{
get
{
return FLoadMomentLeft;
}
set
{
FLoadMomentLeft = value;
}
}
// propertys for rotated data
public double Angle
{
get
{
return FAngle;
}
set
{
FAngle = value;
}
}
public double RemoldFactor
{
get
{
return FRemoldFactor;
}
set
{
FRemoldFactor = value;
}
}
public double SchemaFactor
{
get
{
return FSchemaFactor;
}
set
{
FSchemaFactor = value;
}
}
public double RotateSafetyFactor
{
get
{
return GetRotatedStabilityFactor();
}
}
public bool PC_Ring
{
get
{
return FPC_Ring;
}
set
{
FPC_Ring = value;
}
}
public int WaterlevelNumber
{
get
{
return FWaterlevelNumber;
}
set
{
FWaterlevelNumber = value;
}
}
public int PlotNumber
{
get
{
return FPlotNumber;
}
set
{
FPlotNumber = value;
}
}
public string ErrorNote
{
get
{
return FErrorNote;
}
set
{
FErrorNote = value;
}
}
///
/// Check if the slide plane
/// intersects with the surface
///
public void CheckOnIntersectionWithSurface()
{
var LXIntersect = new List();
var LZIntersect = new List();
double LZ;
FindSurfaceIntersectionPoints(XActiveCentre, ZActiveCentre, ActiveRadius, LXIntersect, LZIntersect);
if (LXIntersect.Count > 1)
{
int i = 0;
// Find entry point (left side )closest to center point
while ((i < LXIntersect.Count))
{
if ((LXIntersect[i] > FxLeftSurface) && (LXIntersect[i] <= XActiveCentre))
{
FErrorNote = FErrorNote + Constants.cIntersectsSurface;
}
i++;
}
}
if ((FErrorNote == ""))
{
LZ = MStabDatafunctions.GetZCoordinateSurfaceAtX(FInterfaceData.SurfaceLine, XActiveCentre);
if ((ZActiveCentre - ActiveRadius) > LZ)
{
FErrorNote = FErrorNote + Constants.cIntersectsSurface;
}
}
if ((FErrorNote == ""))
{
LZ = MStabDatafunctions.GetZCoordinateSurfaceAtX(FInterfaceData.SurfaceLine, XPassiveCentre);
if ((ZPassiveCentre - PassiveRadius) > LZ)
{
FErrorNote = FErrorNote + Constants.cIntersectsSurface;
}
}
if ((FErrorNote == ""))
{
FindSurfaceIntersectionPoints(XActiveCentre, (ZActiveCentre - ActiveRadius), XPassiveCentre,
(ZPassiveCentre - PassiveRadius), LXIntersect, LZIntersect);
}
if ((FErrorNote == ""))
{
LXIntersect.Clear();
LZIntersect.Clear();
FindSurfaceIntersectionPoints(XPassiveCentre, ZPassiveCentre, PassiveRadius, LXIntersect, LZIntersect);
if (LXIntersect.Count > 1)
{
int i = 0;
// Find entry point (left side )closest to center point
while ((i < LXIntersect.Count))
{
if ((LXIntersect[i] < FXRightSurface) && (LXIntersect[i] >= XPassiveCentre))
{
FErrorNote = FErrorNote + Constants.cIntersectsSurface;
}
i++;
}
}
}
}
public override void CreateSlices()
{
TSliceDivision LSliceDivision;
int i;
// Find intersection points with surface
LSliceDivision = new TSliceDivision();
LSliceDivision.CalculationType = FInterfaceData.ModelData.CalcType;
LSliceDivision.LeftSideIsActive = FInterfaceData.ModelData.LeftSideIsActive;
LSliceDivision.StriveSlice = Math.Max(1, FInterfaceData.CalculationOptions.StriveSlice);
LSliceDivision.XLeftSurface = FxLeftSurface;
LSliceDivision.XRightSurface = FXRightSurface;
LSliceDivision.SurfaceLine = FInterfaceData.SurfaceLine;
LSliceDivision.GeometryPoints = FInterfaceData.GeometryPoints;
LSliceDivision.Layers = FInterfaceData.Layers;
LSliceDivision.WaterLines = FInterfaceData.CurrentWaternet.Waterlines;
LSliceDivision.PLLines = FInterfaceData.CurrentWaternet.PLLines;
if (FInterfaceData.CurrentWaternet.PhreaticLineNumber >= 0)
{
LSliceDivision.FreaticLine = FInterfaceData.CurrentWaternet.PLLines[FInterfaceData.CurrentWaternet.PhreaticLineNumber];
}
LSliceDivision.XMidActive = FXMidLeft;
LSliceDivision.ZMidActive = FZMidLeft;
LSliceDivision.XMidpassive = FXMidRight;
LSliceDivision.ZMidpassive = FZMidRight;
LSliceDivision.ZTangent = FZTangent;
LSliceDivision.DivideFailureSurfaceInSlices();
FNSlices = LSliceDivision.NSlices;
FSlices = new Tslice[FNSlices];
for (i = 0; i < FNSlices; i++)
{
FSlices[i] = new Tslice(LSliceDivision.Slices[i]);
}
}
public override double GetReliabilityIndex()
{
double result;
// tijdelijk zolang nog niet in calculation options
bool LIsProbDump;
// LIsReverseGeom: Boolean;
base.GetReliabilityIndex();
result = Constants.CFMinStart;
LIsProbDump = false;
PrepareCalculation();
if ((IsSlipPlaneValid()))
{
// Initialise
CreateSlices();
// find slice numbers from start and end of staight line
FindStraightLineSliceNumbers();
// Fill soil,stress data of slices ( in een apparte routine STSliceFill
FillSlices();
DetermineDrivingMoments();
// nu is bekend welke lagen in het glijvlak zitten, tevens is het goede
// waterspanningsbeeld in de geometrie opgebouwd, Nu met probalistisch
// rekenen (probab of form) beginen
// eerst de variabelen vullen en bijhouden
FStabProbabilisticData = new TMStabProbabilisticData();
try
{
FStabProbabilisticData.InterfaceData = FInterfaceData;
FStabProbabilisticData.Slices = FSlices;
FStabProbabilisticData.CreateandFillProbArrays();
// Shortcuts
FNStochastic = FStabProbabilisticData.NStochastic;
FWithCorrelation = FStabProbabilisticData.WithCorrelation;
FVerdeling = FStabProbabilisticData.Verdeling;
FEx = FStabProbabilisticData.Ex;
FSx = FStabProbabilisticData.Sx;
FVp = FStabProbabilisticData.Vp;
FParamName = FStabProbabilisticData.AParamName;
FCovariantie = FStabProbabilisticData.Covariantie;
FEivec = FStabProbabilisticData.Eivec;
Fx = FStabProbabilisticData.X;
FU = FStabProbabilisticData.U;
FVarArr = FStabProbabilisticData.VarArr;
FAanpasArr = FStabProbabilisticData.AanpasArr;
if (FWithCorrelation)
{
FNumCor = FNStochastic;
FNumEivec = FNStochastic;
}
else
{
FNumCor = 1;
FNumEivec = 1;
}
FProbCalculation = new TProbCalculation();
FProbCalculation.IsTestDump = LIsProbDump;
FProbCalculation.TestFileName = FTestFileName;
FProbCalculation.ZFunction = ZFunction;
FProbCalculation.IsMonteCarlo = FIsMonteCarlo;
FProbCalculation.IsForm = FIsForm;
FProbCalculation.IsDars = FIsDars;
// ResetGlobals(Ex, FVarArr, FAanpasArr);
if ((FIerr == 0))
{
// FStabProbabilisticData.ProbabilisticResultDump(FNStochastic, FBeta, LLevel, LLevelName,
// FX, FU, FSX, FVarArr, FAanpasArr);
// Write output for PC-Ring if appropriate
// if (FPC_Ring) then
// PC_RingResultOutput(FPlotNumber, FWaterlevelNumber, FNStochastic, FBeta, FU, FSX,
// FVarArr, FAanpasArr);
}
// (*
//
// {get the apropriate water level}
// if (FWaterlevelNumber = 0) then
// begin
// LLevel := MyGlobals.StochasticWaterList.MHW
// end
// else
// begin
// LLevel := MyGlobals.StochasticWaterList.stochasticWater[FWaterlevelNumber -
// 1].Level;
// end;
// ShortOutput(Beta, LLevel);
//
// if (not GSTCalculationSettings.ShortReportOn) then
// begin
// WriteMeanAndDesign(Ex, Sx, X, Verdeling, LParamName);
//
// LongOutputProbabilistic(NTot, Beta, LLevel, LLevelName,
// X, U, SX, FVarArr, FAanpasArr);
// end;
//
// Result := Beta;
//
// if (FWaterlevelNumber = 0) then
// begin
// LLevelName := 'Design level';
// end
// else
// begin
// LLevelName := MyGlobals.StochasticWaterList.
// stochasticWater[FWaterlevelNumber - 1].ItemName;
// end;
// {Write results from probabilistic calculation to dump}
// ProbabilisticResultDump(NTot, Beta, LLevel, LLevelName,
// X, U, SX, FVarArr, FAanpasArr);
//
// {Write output for PC-Ring if appropriate}
// if (FPC_Ring) then
// PC_RingResultOutput(FPlotNumber, FWaterlevelNumber, NTot, Beta, U, SX,
// FVarArr, FAanpasArr);
// end;
// *)
// Free dynamic arrays
}
finally
{
FVarArr = null;
FAanpasArr = null;
FVerdeling = null;
FEx = null;
FSx = null;
FVp = null;
FCovariantie = null;
FEivec = null;
Fx = null;
FU = null;
FParamName = null;
}
}
return result;
}
protected void FindSurfaceIntersectionPoints(double AX1, double AZ1, double AX2, double AZ2, List AXIntersect, List AZIntersect)
{
// find intersection points
for (int i = FInterfaceData.SurfaceLine.GetLowerBound(0) + 1; i <= FInterfaceData.SurfaceLine.GetUpperBound(0); i++)
{
double XIntersect = 0;
double ZIntersect = 0;
double LXStart = FInterfaceData.SurfaceLine[i - 1].XCoor;
double LZStart = FInterfaceData.SurfaceLine[i - 1].ZCoor;
double LXEnd = FInterfaceData.SurfaceLine[i].XCoor;
double LZEnd = FInterfaceData.SurfaceLine[i].ZCoor;
if (MStabDatafunctions.IntersectLines(AX1, AZ1, AX2, AZ2, LXStart, LZStart, LXEnd, LZEnd, ref XIntersect, ref ZIntersect))
{
FErrorNote = FErrorNote + Constants.cIntersectsSurface;
}
}
}
// =======================================================================================================================
// Date ID Modification
// 2008-06-25 Best Created
// =======================================================================================================================
protected void FindSurfaceIntersectionPoints(double AXMid, double AZMid, double ARad, List AXIntersect, List AZIntersect)
{
// find intersection points
for (int i = FInterfaceData.SurfaceLine.GetLowerBound(0) + 1; i <= FInterfaceData.SurfaceLine.GetUpperBound(0); i++)
{
double LXStart = FInterfaceData.SurfaceLine[i - 1].XCoor;
double LZStart = FInterfaceData.SurfaceLine[i - 1].ZCoor;
double LXEnd = FInterfaceData.SurfaceLine[i].XCoor;
double LZEnd = FInterfaceData.SurfaceLine[i].ZCoor;
MStabDatafunctions.SnijCirkelLijnstuk(AXMid, AZMid, ARad, LXStart, LZStart, LXEnd, LZEnd, AXIntersect, AZIntersect);
}
MStabDatafunctions.SortIntersectionPoints(AZIntersect, AXIntersect);
MStabDatafunctions.RemoveDoubleIntersectionPoints(AXIntersect, AZIntersect);
}
// =======================================================================================================================
// Date ID Modification
// 2008-06-23 Best Created
// =======================================================================================================================
protected void FindEntreepoints()
{
var LXIntersect = new List();
var LZIntersect = new List();
bool LPointTooHigh = false;
bool LPointLeft = false;
FxLeftSurface = Double.NaN;
FZLeftSurface = Double.NaN;
FindSurfaceIntersectionPoints(FXMidLeft, FZMidLeft, FRadLeft, LXIntersect, LZIntersect);
if (LXIntersect.Count == 1 && LXIntersect[0] <= FXMidLeft)
{
FxLeftSurface = LXIntersect[0];
FZLeftSurface = LZIntersect[0];
LPointTooHigh = (FZLeftSurface > FZMidLeft);
}
else if (LXIntersect.Count >= 2)
{
int i = 0;
// Find entry point (left side )closest to center point
while ((i < LXIntersect.Count - 1))
{
if (LXIntersect[i] <= FXMidLeft)
{
LPointLeft = true;
FxLeftSurface = LXIntersect[i];
FZLeftSurface = LZIntersect[i];
LPointTooHigh = (FZLeftSurface > FZMidLeft);
}
i++;
}
}
else
{
FErrorNote = FErrorNote + Constants.cNoEntryPointFound;
}
if (LPointTooHigh)
{
FErrorNote = FErrorNote + Constants.cEntryPointTooHigh;
}
if (!LPointLeft)
{
FErrorNote = FErrorNote + Constants.cNoEntryPointFound;
}
}
// =======================================================================================================================
// Date ID Modification
// 2008-06-23 Best Created
// =======================================================================================================================
protected void FindExitpoints()
{
var LXIntersect = new List();
var LZIntersect = new List();
bool LPointTooHigh = false;
bool LPointRight = false;
FindSurfaceIntersectionPoints(FXMidRight, FZMidRight, FRadRight, LXIntersect, LZIntersect);
FXRightSurface = Double.NaN;
FZRightSurface = Double.NaN;
if (LXIntersect.Count == 1 && LXIntersect[0] >= FXMidRight)
{
FXRightSurface = LXIntersect[0];
FZRightSurface = LZIntersect[0];
LPointTooHigh = (LZIntersect[0] > FZMidRight);
}
else if (LXIntersect.Count >= 2)
{
int i = 1;
bool LFound = false;
// find exit point (Right side) closest to center point
while ((i <= LXIntersect.Count - 1) && !LFound)
{
LFound = (LXIntersect[i] >= FXMidRight);
if (LFound)
{
LPointRight = true;
FXRightSurface = LXIntersect[i];
FZRightSurface = LZIntersect[i];
LPointTooHigh = (LZIntersect[i] > FZMidRight);
}
i++;
}
}
else
{
FErrorNote = FErrorNote + Constants.cNoExitPointFound;
}
if (LPointTooHigh)
{
FErrorNote = FErrorNote + Constants.cExitPointTooHigh;
}
if (!LPointRight)
{
FErrorNote = FErrorNote + Constants.cNoExitPointFound;
}
}
// =======================================================================================================================
// Date ID Modification
// 2008-06-25 Best Created
// =======================================================================================================================
protected void CheckOnForbiddenLines()
{
int i;
double LRadius;
int LNDim;
bool LIntersects;
LIntersects = false;
LNDim = FInterfaceData.ForbiddenLines.Length;
if ((LNDim > 0))
{
i = 0;
while ((i < LNDim) && !LIntersects)
{
TForbiddenLine _wvar1 = FInterfaceData.ForbiddenLines[i];
// halve cirkel aan linkerkant
LRadius = FZMidLeft - FZTangent;
if (!LIntersects)
{
LIntersects = _wvar1.CircleIntersectsForbiddenLine(FXMidLeft, FZMidLeft, LRadius, FxLeftSurface, FXMidLeft);
}
if (!LIntersects)
{
LIntersects = _wvar1.LineIntersectsForbiddenLine(FXMidLeft, FZTangent, FXMidRight, FZTangent);
}
LRadius = FZMidRight - FZTangent;
if (!LIntersects)
{
LIntersects = _wvar1.CircleIntersectsForbiddenLine(FXMidRight, FZMidRight, LRadius, FXMidRight, FXRightSurface);
}
i++;
}
}
if (LIntersects)
{
FErrorNote = FErrorNote + Constants.cIntersectsForbiddenLine;
}
}
protected bool IsSlipPlaneValid()
{
bool result;
string lStr;
try
{
FindEntreepoints();
FindExitpoints();
if (Double.IsNaN(FxLeftSurface) || Double.IsNaN(FZLeftSurface) ||
Double.IsNaN(FXRightSurface) || Double.IsNaN(FZRightSurface))
{
return false;
}
// See if surface intersects slide plane
if ((FErrorNote == ""))
{
CheckOnIntersectionWithSurface();
}
// See if intersection with forbidden lines are found
if ((FErrorNote == ""))
{
CheckOnForbiddenLines();
}
if ((FErrorNote == ""))
{
FErrorNote = TestCenterpointInGeometry(FXPassiveCentre, FZPassiveCentre);
}
if ((FErrorNote == ""))
{
FErrorNote = TestCenterpointInGeometry(FXActiveCentre, FZActiveCentre);
}
result = (FErrorNote == "");
}
catch (Exception E)
{
lStr = E.Message;
throw E;
}
return result;
}
// internal data
// =======================================================================================================================
// Date ID Modification
// 2008-08-06 Best Created
// =======================================================================================================================
private void PrepareCalculation()
{
FGamWat = FInterfaceData.CurrentWaternet.GammaWater;
if ((FXActiveCentre < FXPassiveCentre))
{
FXMidLeft = FXActiveCentre;
FZMidLeft = FZActiveCentre;
FRadLeft = FZActiveCentre - FZTangent;
FXMidRight = FXPassiveCentre;
FZMidRight = FZPassiveCentre;
FRadRight = FZPassiveCentre - FZTangent;
}
else
{
FXMidLeft = FXPassiveCentre;
FZMidLeft = FZPassiveCentre;
FRadLeft = FZPassiveCentre - FZTangent;
FXMidRight = FXActiveCentre;
FZMidRight = FZActiveCentre;
FRadRight = FZActiveCentre - FZTangent;
}
}
// function GetYAtCircle(thex, xmid, ymid, radius: Double): Double;
// =======================================================================================================================
// Description: Calculates the DrivingSoilMoment on the left side
//
// Date ID Modification
// 1999-03-29 Best Created
// =======================================================================================================================
private double DrivingSoilMomentLeft()
{
double result;
double HorDistance;
double Moment;
int i;
Moment = 0;
for (i = 0; i <= FLeftNumber; i++)
{
Tslice _wvar1 = FSlices[i];
HorDistance = FXMidLeft - _wvar1.XMid;
// Bereken aandrijvend moment : dx * (lamel gewicht - water gewicht)
Moment = Moment + HorDistance*_wvar1.Weight;
}
result = Moment;
return result;
}
// =======================================================================================================================
// Description: Calculates the DrivingSoilMoment on the right side
//
// Date ID Modification
// 1999-03-29 Best Created
// =======================================================================================================================
private double DrivingSoilMomentRight()
{
double result;
double HorDistance;
double Moment;
int i;
Moment = 0;
for (i = FRightNumber; i < FNSlices; i++)
{
Tslice _wvar1 = FSlices[i];
HorDistance = FXMidRight - _wvar1.XMid;
// Bereken aandrijvend moment : dx * (lamel gewicht - water gewicht)
Moment = Moment + HorDistance*_wvar1.Weight;
}
result = Moment;
return result;
}
// =======================================================================================================================
// Date ID Modification
// 2008-08-07 Best Created
// =======================================================================================================================
private double WaterMomentOnVerticals(double AXleftside, double AXRightside, double AXCentre, double AZCentre, double ARadius)
{
double result = 0;
double Ztop = 0;
double ZBot = 0;
double Ptop = 0;
double PBot = 0;
double ZFrea = 0;
double dZ = 0;
double force = 0;
double LMom = 0;
double LWaterMoment = 0;
bool LDown = false;
double LZFreaTop = 0;
double LZFreaBottom = 0;
int LNS = 0;
double LXs1 = 0;
double LXS2 = 0;
double Lzs1 = 0;
double LzS2 = 0;
bool LVerticalFound = false;
double LX = 0;
double LZ1 = 0;
double LZ2 = 0;
TPoints[] LFreaticLine;
if (FInterfaceData.CurrentWaternet.PhreaticLineNumber == -1)
{
return 0;
}
LFreaticLine = FInterfaceData.CurrentWaternet.PLLines[FInterfaceData.CurrentWaternet.PhreaticLineNumber];
LWaterMoment = 0;
MStabDatafunctions.FindVertical(FInterfaceData.SurfaceLine, ref LVerticalFound, ref LX, ref LZ1, ref LZ2);
if ((LVerticalFound && (LX >= AXleftside - Constants.CSameSliceDist) && (LX <= AXRightside + Constants.CSameSliceDist) && (LZ1 != LZ2)))
{
LDown = (LZ1 > LZ2);
// calculate freatop and bottom in case of a sheetpiling
LZFreaTop = MStabDatafunctions.ZTopAtSurface(LX, LFreaticLine);
LZFreaBottom = MStabDatafunctions.ZBottomAtSurface(LX, LFreaticLine);
if ((Math.Abs(LZFreaTop - LZFreaBottom) < Constants.CGeoAccu))
{
ZFrea = LZFreaTop;
}
else
{
// Find the correct freatic line in case of a jump in freatic lines
if (LDown)
{
ZFrea = MStabDatafunctions.ZTopAtSurface(LX + Constants.CGeoAccu, LFreaticLine);
}
else
{
ZFrea = MStabDatafunctions.ZTopAtSurface(LX - Constants.CGeoAccu, LFreaticLine);
}
}
Ztop = Math.Max(LZ1, LZ2);
ZBot = Math.Min(LZ1, LZ2);
// see if vertical line is in reach of circle
LMom = 0;
if ((MStabDatafunctions.Distance2D(LX, ZBot, AXCentre, AZCentre) < ARadius) || (MStabDatafunctions.Distance2D(LX, Ztop, AXCentre, AZCentre) < ARadius))
{
// Check if entry/exit point is on the curve
if ((MStabDatafunctions.Distance2D(LX, ZBot, AXCentre, AZCentre) > ARadius))
{
MStabDatafunctions.Intersect_Circle_line(AXCentre, AZCentre, ARadius, LX, LX, ZBot, Ztop, ref LNS, ref LXs1, ref Lzs1, ref LXS2, ref LzS2);
if (LNS == 1)
{
if (Lzs1 >= ZBot && Lzs1 <= Ztop)
{
ZBot = Lzs1;
}
}
else if ((LNS == 2))
{
if ((Lzs1 >= ZBot) && (Lzs1 <= Ztop))
{
ZBot = Lzs1;
}
else
{
ZBot = LzS2;
}
}
}
dZ = Ztop - ZBot;
// Case 1: Water level above top point of vertical ->
// Calculate trapezium shaped water pressure figure
if (ZFrea > Ztop)
{
Ptop = (ZFrea - Ztop)*FGamWat;
PBot = (ZFrea - ZBot)*FGamWat;
force = Ptop*(Ztop - ZBot);
// Rectangular part
LMom = force*(AZCentre - (ZBot + dZ/2));
force = 0.5*dZ*(PBot - Ptop);
// Triangular part
LMom = LMom + force*(AZCentre - (ZBot + dZ/3));
}
else
{
// Case 2: Water level intersects vertical surface piece.
// Calculate triangular shaped water pressure figure
if (ZFrea > ZBot)
{
PBot = (ZFrea - ZBot)*FGamWat;
force = 0.5*(ZFrea - ZBot)*PBot;
LMom = force*(AZCentre - (ZBot + (ZFrea - ZBot)/3));
}
else
{
// Case 3: Water level below vertical surface piece -> M=0
LMom = 0;
}
}
// If vertical line goes 'down', the water force points
// 'to the left' therefore the moment gets a '-' sign
if (LDown)
{
LMom = -LMom;
}
}
// Sum contributions
LWaterMoment = LWaterMoment + LMom;
}
// Return total moment
result = LWaterMoment;
return result;
}
// =======================================================================================================================
// Description: A uplift Van calculation consist of two circel parts and a horicontal part
// Equilibrium for the circel parts must include Water against such a part
//
// Date ID Modification
// 2007-02-15 Best Created
// =======================================================================================================================
private double WaterMomentOnCircleEnd(double AXCoor, double AXCentre, double AZCentre, double ARadius, bool AIsLeft)
{
double result;
double LWaterMoment;
double LPoretop;
double LPoreBot;
double LRectPore;
double LForce;
double LDz;
double LXCoor;
// LVerticalData: TVerticalData;
double LZCoor;
double LLastZCoor;
double LZ = 0;
double LZHead = 0;
TVertical LVertical;
double LSoilTop;
double LFreaticLevel;
LZCoor = AZCentre - ARadius;
LXCoor = AXCoor;
LWaterMoment = 0;
LVertical = new TVertical();
LVertical.XCoordinate = LXCoor;
LVertical.ZCoordinate = LZCoor;
LVertical.SoilData.AddRange(FSoilData);
LVertical.GammaWater = FInterfaceData.CurrentWaternet.GammaWater;
if (FInterfaceData.CurrentWaternet.PhreaticLineNumber != -1)
{
LVertical.FreaticLine = FInterfaceData.CurrentWaternet.PLLines[FInterfaceData.CurrentWaternet.PhreaticLineNumber];
}
LVertical.WaterLines.AddRange(FInterfaceData.CurrentWaternet.Waterlines);
LVertical.PLLines = FInterfaceData.CurrentWaternet.PLLines;
LVertical.FillVerticalWithWaterData();
LVertical.SurfaceLine = FInterfaceData.SurfaceLine;
LVertical.InterfaceCalls = FInterfaceData.InterfaceCalls;
LVertical.FillVertical();
LSoilTop = LVertical.SoilLevel;
LFreaticLevel = LVertical.Freaticlevel;
LPoreBot = LVertical.DetermineWaternetStress(LZCoor);
LLastZCoor = LZCoor;
// find net waterline head below surface
while (LVertical.NextWaterLineIsFound(LLastZCoor, ref LZ, ref LZHead))
{
LPoretop = (LZHead - LZ)*FGamWat;
// rectangualr moment if present
LRectPore = Math.Min(LPoretop, LPoreBot);
LDz = LZ - LLastZCoor;
LForce = LRectPore*LDz;
LWaterMoment = LWaterMoment + LForce*(AZCentre - (LLastZCoor + LDz/2));
// triangulair part if present
LForce = 0.5*LDz*Math.Abs(LPoretop - LPoreBot);
if ((LPoretop < LPoreBot))
{
LDz = LLastZCoor + LDz/3;
}
else
{
LDz = LZ - LDz/3;
}
LWaterMoment = LWaterMoment + LForce*(AZCentre - LDz);
LLastZCoor = LZ;
LPoreBot = LPoretop;
}
// Remove moment above surface as this is already teken into acount
if ((LFreaticLevel > LSoilTop))
{
LDz = LFreaticLevel - LSoilTop;
LForce = 0.5*LDz*LDz*FInterfaceData.CurrentWaternet.GammaWater;
LWaterMoment = LWaterMoment - LForce*(AZCentre - LSoilTop + LDz/3);
}
if (AIsLeft)
{
result = -LWaterMoment;
}
else
{
result = LWaterMoment;
}
return result;
}
// =======================================================================================================================
// Description: Calculates the DrivingWaterMoment
//
// Date ID Modification
// 1997-07-01 Best Created
// =======================================================================================================================
private double DrivingWaterMomentLeft()
{
double result;
double HorDistance;
double VertDistance;
double LWaterMom;
double Moment;
int i;
bool LIsLeft = false;
Moment = 0;
for (i = 0; i <= FLeftNumber; i++)
{
Tslice _wvar1 = FSlices[i];
// Calculate watermoment due to waterforces on slices
HorDistance = _wvar1.XMid - FXMidLeft;
VertDistance = FZMidLeft - _wvar1.ZTopMid;
Moment = Moment + HorDistance*FSlices[i].VPoreOnSurface + VertDistance*FSlices[i].HPoreOnSurface;
}
// Calculate extra moments due to water forces on vertical parts of
// the geometry within the slip circle area.
Moment = Moment + WaterMomentOnVerticals(FxLeftSurface, FXMidLeft, FXMidLeft, FZMidLeft, FZMidLeft - FZTangent);
// Calculate extra moments due to water forces on vertical parts of
// the geometry within the slip circle area.
LIsLeft = true;
LWaterMom = WaterMomentOnCircleEnd(FXMidLeft, FXMidLeft, FZMidLeft, FZMidLeft - FZTangent, LIsLeft);
Moment = Moment + LWaterMom;
result = Moment;
return result;
}
// --- Einde procedure Momenten ---
// =======================================================================================================================
// Method : DrivingWaterMoment : Double;
//
// Member of class : TSTLiftPlane.
//
// Acces : Public
//
// Description: Calculates the DrivingWaterMoment
//
// Pre-condition : The Circle,Geometry must exist and slices are filled.
//
// Post-condition :
//
// Method History :
//
// Date ID Modification
// 1997-07-01 Best Created
// =======================================================================================================================
private double DrivingWaterMomentRight()
{
double result;
double HorDistance;
double VertDistance;
double Moment;
int i;
bool LIsLeft;
Moment = 0;
for (i = FRightNumber; i < FNSlices; i++)
{
Tslice _wvar1 = FSlices[i];
// Calculate watermoment due to waterforces on slices
HorDistance = _wvar1.XMid - FXMidRight;
VertDistance = FZMidRight - _wvar1.ZTopMid;
Moment = Moment + HorDistance*_wvar1.VPoreOnSurface + VertDistance*_wvar1.HPoreOnSurface;
}
// Calculate extra moments due to water forces on vertical parts of
// the geometry within the slip circle area.
Moment = Moment + WaterMomentOnVerticals(FXMidRight, FXRightSurface, FXMidRight, FZMidRight, FZMidRight - FZTangent);
LIsLeft = false;
Moment = Moment + WaterMomentOnCircleEnd(FXMidRight, FXMidRight, FZMidRight, FZMidRight - FZTangent, LIsLeft);
result = Moment;
return result;
}
// --- Einde procedure Momenten ---
// =======================================================================================================================
// Description: Calculates the moments of the non uniform loads
//
// Date ID Modification
// 1999-03-09 Best Created
// =======================================================================================================================
private double UniformLoadMoment(double AXCentre, double AXLeft, double AXRight)
{
double result;
int i;
double Som;
Som = 0.0;
for (i = FInterfaceData.UniformLoad.GetLowerBound(0); i <= FInterfaceData.UniformLoad.GetUpperBound(0); i++)
{
TUniformLoad _wvar1 = FInterfaceData.UniformLoad[i];
Som = Som + _wvar1.GetMomentRoundX(AXCentre, AXLeft, AXRight);
}
result = Som;
return result;
}
// =======================================================================================================================
// Date ID Modification
// 2008-08-07 Best Created
// =======================================================================================================================
private double LineLoadMoment(double AXCentre, double AZCentre, double ARadius, double AXLeft, double AXRight)
{
double result;
int i;
double Som;
Som = 0.0;
for (i = FInterfaceData.LineLoad.GetLowerBound(0); i <= FInterfaceData.LineLoad.GetUpperBound(0); i++)
{
TLineLoad _wvar1 = FInterfaceData.LineLoad[i];
Som = Som + _wvar1.GetMoment(AXCentre, AZCentre, ARadius, AXLeft, AXRight);
}
result = Som;
return result;
}
// =======================================================================================================================
// Date ID Modification
// 2008-08-06 Best Created
// =======================================================================================================================
private void DetermineDrivingMoments()
{
FTotDrivingLeft = 0.0;
FTotDrivingRight = 0.0;
// calculate driving moments
FSoilMomentLeft = DrivingSoilMomentLeft();
FSoilMomentRight = DrivingSoilMomentRight();
// watermoments from free water
FWatermomentLeft = DrivingWaterMomentLeft();
FWatermomentRight = DrivingWaterMomentRight();
// load moments if present
FLoadMomentLeft = UniformLoadMoment(FXMidLeft, FxLeftSurface, FXMidLeft) + LineLoadMoment(FXMidLeft, FZMidLeft, FRadLeft, FxLeftSurface, FXMidLeft);
FLoadMomentRight = UniformLoadMoment(FXMidRight, FXMidRight, FXRightSurface) + LineLoadMoment(FXMidRight, FZMidRight, FRadRight, FXMidRight, FXRightSurface);
// the total driving moments
FTotDrivingLeft = FSoilMomentLeft + FWatermomentLeft + FLoadMomentLeft;
FTotDrivingRight = FSoilMomentRight + FWatermomentRight + FLoadMomentRight;
FReverseGeom = Math.Abs(FTotDrivingRight) > Math.Abs(FTotDrivingLeft);
}
// =======================================================================================================================
// Description: Calculates the horizontal force between two circle parts in MLift
//
// Date ID Modification
// 1999-03-10 Best Created
// 2007-03-05 Best rule implemented for horizontal part too
// =======================================================================================================================
private double GetHorResistingForce(int StartLam, int EndLam, double Fo)
{
double result = 0;
const double Pir = Math.PI/180.0;
int i = 0;
int TsPoint = 0;
double Som = 0;
double coh = 0;
double phi = 0;
double tgPhi = 0;
double LCohTerm = 0;
double LSigTerm = 0;
double LSinPhi = 0;
double LCosPhi = 0;
double LCosDilet = 0;
double LSinDilet = 0;
double LTau = 0;
Som = 0.0;
for (i = StartLam; i <= EndLam; i++)
{
Tslice _wvar1 = FSlices[i];
coh = _wvar1.CohBottom;
// cohesie
LSinPhi = Math.Sin(GetSlice(i).PhiBottom*Pir);
LCosPhi = Math.Cos(_wvar1.PhiBottom*Pir);
LCosDilet = Math.Cos(_wvar1.Dilatancy*Pir);
LSinDilet = Math.Sin(_wvar1.Dilatancy*Pir);
LSigTerm = LCosDilet*LSinPhi/(1 - LSinDilet*LSinPhi);
LCohTerm = LCosDilet*LCosPhi/(1 - LSinDilet*LSinPhi);
// bereken de coh en phi uit de spannings tabel
if (new ArrayList(new TMatStrengthTypeSet[]
{
TMatStrengthTypeSet.mstStressTab, TMatStrengthTypeSet.mstPseudoStressTab
}).Contains(FSlices[i].StrengthType))
{
TStressTable _wvar2 = InterfaceData.Stresstable[_wvar1.SpanTabel];
_wvar2.GetCohPhi(_wvar1.EffectiveStress, _wvar1.PseudoFactor, ref coh, ref phi, ref TsPoint);
_wvar1.PhiBottom = phi;
// New phi - alfa
_wvar1.CohBottom = coh;
// New coh - alfa
tgPhi = Math.Tan(_wvar1.PhiBottom*Pir);
// New tangens phie
LCohTerm = 1;
LSigTerm = tgPhi;
}
LTau = coh*LCohTerm + FSlices[i].EffectiveStress*LSigTerm;
Som = Som + _wvar1.ArcLength*LTau;
// Nieuwe Tau - alfa
_wvar1.NormalStress = _wvar1.EffectiveStress;
FSlices[i].ShearStress = LTau/Fo;
}
result = Som;
return result;
}
// =======================================================================================================================
// Description: Calculates the resisting moments of a circle part
//
// Date ID Modification
// 1999-03-31 Best Created
// =======================================================================================================================
private double GetResistingMoment(int Startlam, int Endlam, double Radius, double fo)
{
double result = 0;
const double Pir = Math.PI/180.0;
int i = 0;
int TsPoint = 0;
double coh = 0;
double phi = 0;
double tgPhi = 0;
double Salfa = 0;
double alfa = 0;
double h = 0;
double tgalfa = 0;
double SomMoment = 0;
double LCohTerm = 0;
double LSigTerm = 0;
double LSinPhi = 0;
double LCosPhi = 0;
double LCosDilet = 0;
double LSinDilet = 0;
double LTau = 0;
SomMoment = 0.0;
for (i = Startlam; i <= Endlam; i++)
{
coh = FSlices[i].CohBottom;
// cohesie
LSinPhi = Math.Sin(FSlices[i].PhiBottom*Pir);
LCosPhi = Math.Cos(FSlices[i].PhiBottom*Pir);
LCosDilet = Math.Cos(FSlices[i].Dilatancy*Pir);
LSinDilet = Math.Sin(FSlices[i].Dilatancy*Pir);
tgPhi = Math.Tan(FSlices[i].PhiBottom*Pir);
// formula based on Tau-alfa =
// 1/Fo * (coh * cotgphi + Salfa')* cospsi*sinphi/(1 - sinphi*sinpsi)
// (after tuenissen))
// note that if psi = phi then tau = coh + salfa * tgphi
LSigTerm = LCosDilet*LSinPhi/(1 - LSinDilet*LSinPhi);
LCohTerm = LCosDilet*LCosPhi/(1 - LSinDilet*LSinPhi);
alfa = FSlices[i].BottomAngle;
// Hoek onderkant lamel [rad]
if (FDrivingMoment > 0.0)
{
alfa = -alfa;
}
// Passieve alfa nagatief maken
// Afsnuiten indien nodig
h = 0.5*Math.Atan(tgPhi/fo) - 0.25*Math.PI;
if (alfa <= h)
{
alfa = h;
}
tgalfa = Math.Tan(alfa);
// Bereken het tegen werkend moment
LTau = coh*LCohTerm + FSlices[i].EffectiveStress*LSigTerm;
SomMoment = SomMoment + FSlices[i].ArcLength*LTau/(1 + tgalfa*LSigTerm/fo);
Salfa = (fo*FSlices[i].EffectiveStress - coh*LCohTerm*tgalfa)/(fo + LSigTerm*tgalfa);
FSlices[i].NormalStress = Salfa;
FSlices[i].ShearStress = (coh*LCohTerm + Salfa*LSigTerm)/fo;
// (*
// { Calculate resisting moment }
// SomMoment := SomMoment + ArcLength * (coh + EffectiveStress * tgphi) /
// (1 + tgalfa * tgphi / Fo);
//
// { Vertical balance : Sigma-v' = Salfa' + Tau-alfa*Tgalfa }
// { en Tau according to Bishop : Tau-alfa = 1/Fo * (coh + Salfa'*tgphi }
// { From these 2 equations Salfa can be solved : }
// Salfa := (Fo * EffectiveStress - coh * tgalfa) / (Fo + tgalfa * tgphi);
//
// NormalStress := Salfa;
// ShearStress := (Coh + Salfa * tgphi) / Fo;
// *)
if (new ArrayList(new TMatStrengthTypeSet[]
{
TMatStrengthTypeSet.mstStressTab, TMatStrengthTypeSet.mstPseudoStressTab
}).Contains(FSlices[i].StrengthType))
{
// Nieuwe sigma alfa
if (Salfa < 0.0)
{
Salfa = 0.0;
// Add message only if not already added earlier
if (ErrorNote.IndexOf("Negative effective stress.") == 0)
{
ErrorNote = ErrorNote + "Negative effective stress.@";
}
}
FSlices[i].SigmaAlfa = Salfa;
FInterfaceData.StabSoilData[FSlices[i].SoilNumber].GetCohAndPhiFromTable(FSlices[i].SigmaAlfa,
FSlices[i].PseudoFactor, ref coh, ref phi, ref TsPoint);
FSlices[i].PhiBottom = phi;
// New phi - alfa
// dilatancy wordt alleen in methode c,phi gebruikt
FSlices[i].Dilatancy = phi;
FSlices[i].CohBottom = coh;
// New coh - alfa
tgPhi = Math.Tan(FSlices[i].PhiBottom*Pir);
// New tangens phie
// Nieuwe Tau - alfa
FSlices[i].NormalStress = Salfa;
FSlices[i].ShearStress = (FSlices[i].CohBottom + Salfa*tgPhi)/fo;
}
}
result = SomMoment*Radius;
return result;
}
// =======================================================================================================================
// Date ID Modification
// 2008-03-26 Best Created
// =======================================================================================================================
// Internal utility methods
// procedure RemoveLocalPointers;
// function WaterMomentOnVerticals(Xleftside, XRightside: Double;
// XCentre, YCentre, Radius: Double): Double;
// Private Declarations
// =======================================================================================================================
// Date ID Modification
// 2008-08-06 Best Created
// =======================================================================================================================
private void FindStraightLineSliceNumbers()
{
double LXMLeft;
double LXMRight;
int i;
// Find the last and first slicenumbers of the cirkels
// het rechte stuk loopt dus van leftnumber + 1 tot rightnumber - 1
FLeftNumber = 0;
FRightNumber = FNSlices - 1;
LXMLeft = Math.Min(FXActiveCentre, FXPassiveCentre);
LXMRight = Math.Max(FXActiveCentre, FXPassiveCentre);
for (i = 0; i < FNSlices; i++)
{
if ((Math.Abs(LXMLeft - FSlices[i].xRight) < Constants.CSameSliceDist))
{
FLeftNumber = i;
}
if ((Math.Abs(LXMRight - FSlices[i].xLeft) < Constants.CSameSliceDist))
{
FRightNumber = i;
}
}
}
// =======================================================================================================================
//
// Procedure : LiftMethode (Iter)
//
// Declared in unit : TSTLiftPlane (Internal)
//
// Description : Berekent de veiligheidscoeff. uit bij methode
// MLift.
//
// Last Update 19990329 : (Best) Created
//
// Name Type Function
// ---- ---- --------
// Parameters - IN : Iter Integer Houdt antal iteraties bij
//
// Pre-condition : None
//
// Post-condition : None
// 2007-02-15 Best
// =======================================================================================================================
private double LiftMethode(ref int AIter)
{
double result;
double ArmLeft;
double ArmRight;
double LeftResMom;
double RightResMom;
double HorResForce;
double Fn;
double Fo;
bool Ready;
double LHorWaterForce;
// Initialisatie
FResistingMom = 0.0;
Fn = FInterfaceData.CalculationOptions.StartValueSafetyFactor;
Ready = false;
AIter = 0;
Tslice _wvar1 = FSlices[FLeftNumber];
ArmLeft = FRadLeft - 1d/3d*(_wvar1.ZTopRight - _wvar1.ZBottomRight);
Tslice _wvar2 = FSlices[FRightNumber];
ArmRight = FRadRight - 1d/3d*(_wvar2.ZTopLeft - _wvar2.ZBottomLeft);
FDrivingMoment = FTotDrivingLeft/ArmLeft + FTotDrivingRight/ArmRight;
// Driving moment is now actualy a force
LHorWaterForce = DetermineHoriZontalWaterforce();
//add the resultant of the water force
FDrivingMoment = FDrivingMoment + LHorWaterForce;
FTextileMoment = GeotextileMoment(FDrivingMoment);
// Check driving moment on 1 kNm
if ((Math.Abs(FDrivingMoment) < 0.001))
{
Ready = true;
}
// Start iteratie
if (!Ready)
{
do
{
AIter++;
Fo = Fn;
FResistingMom = 0.0;
// bereken de momenten van de linker en rechter helfr
// rekening houdend met de draairichting (reversegeom of niet )
LeftResMom = GetResistingMoment(0, FLeftNumber, FRadLeft, Fo);
RightResMom = GetResistingMoment(FRightNumber, FNSlices - 1, FRadRight, Fo);
HorResForce = GetHorResistingForce(FLeftNumber + 1, FRightNumber - 1, Fo);
FResistingMom = Math.Abs(LeftResMom)/ArmLeft + Math.Abs(RightResMom)/ArmRight + Math.Abs(HorResForce);
// bereken de actieve passieve en hor kracht
if (FReverseGeom)
{
FPassiveForce = (Math.Abs(FTotDrivingLeft) + Math.Abs(LeftResMom/Fo))/ArmLeft;
FActiveForce = (Math.Abs(FTotDrivingRight) - Math.Abs(RightResMom/Fo))/ArmRight;
FActiveDMoment = FTotDrivingRight;
FActiveRMoment = RightResMom;
FPassiveDMoment = FTotDrivingLeft;
FPassiveRMoment = LeftResMom;
}
else
{
FActiveForce = (Math.Abs(FTotDrivingLeft) - Math.Abs(LeftResMom)/Fo)/ArmLeft;
FPassiveForce = (Math.Abs(FTotDrivingRight) + Math.Abs(RightResMom/Fo))/ArmRight;
FActiveDMoment = FTotDrivingLeft;
FActiveRMoment = LeftResMom;
FPassiveDMoment = FTotDrivingRight;
FPassiveRMoment = RightResMom;
}
FHorForce = HorResForce/Fo;
if ((AIter == 1))
{
FActiveForce0 = FActiveForce;
FPassiveForce0 = FPassiveForce;
FHorForce0 = FHorForce;
}
FResistingMom = FResistingMom + FTextileMoment;
Fn = FResistingMom/Math.Abs(FDrivingMoment);
// Sla ongeitereerd moment, dus bij 1e iteratie op, om te kunnen uitvoeren
if (AIter == 1)
{
FResistingMom0 = FResistingMom;
}
// Stoppen als veiligheidscoefficient oud en nieuw ongeveer gelijk zijn
if (Math.Abs(Fn - Fo) < 0.001)
{
Ready = true;
}
// Teveel iteratie dan stoppen
if (AIter > Constants.CMaxIter)
{
Ready = true;
}
// Als tegenwerkend moment = 0, dan ook stoppen.
if ((Math.Abs(FResistingMom) < 0.001) || (Fn < 0.001))
{
Ready = true;
}
} while (!(Ready));
}
result = Fn;
if (Math.Abs(FDrivingMoment) < 0.001)
{
result = 1001;
LogManager.Add(new LogMessage(LogMessageType.Error, null, LocalizationManager.GetTranslatedText(GetType(), "DrivingMomentTooSmall")));
ErrorNote = ErrorNote + "Driving moment too small.@";
}
else
{
if (AIter >= Constants.CMaxIter)
{
result = 1001;
LogManager.Add(new LogMessage(LogMessageType.Error, null, string.Format(LocalizationManager.GetTranslatedText(GetType(), "TooManyIterations"), Constants.CMaxIter)));
FErrorNote = FErrorNote + "More than " + Convert.ToString(Constants.CMaxIter) + " iterations.@";
}
else if ((Math.Abs(FResistingMom) < 0.001) || (Fn < 0.001))
{
result = 1001;
LogManager.Add(new LogMessage(LogMessageType.Error, null, LocalizationManager.GetTranslatedText(GetType(), "ResistingMomentTooSmall")));
FErrorNote = FErrorNote + "Resisting moment too small.@";
}
}
return result;
}
private double determinePoreOnSurface(double aXCoor)
{
double result = 0;
double freaticlevel;
if (FInterfaceData.CurrentWaternet.PhreaticLineNumber != -1)
{
TPoints[] FreaticLine =
FInterfaceData.CurrentWaternet.PLLines[FInterfaceData.CurrentWaternet.PhreaticLineNumber];
freaticlevel = MStabDatafunctions.ZTopAtSurface(aXCoor, FreaticLine);
TPoints[] SurfaceLine = FInterfaceData.SurfaceLine;
double lTop = MStabDatafunctions.ZTopAtSurface(aXCoor, SurfaceLine);
if (freaticlevel > lTop)
{
result = (freaticlevel - lTop)*FInterfaceData.CurrentWaternet.GammaWater;
result = result*0.5*(freaticlevel - lTop);
}
}
return result;
}
private double DetermineHoriZontalWaterforce()
{
double LLeftForce;
double LRightForce;
double LRadius;
double LeftPoreOnSurface;
double RightPoreOnSurface;
LRadius = FZMidLeft - FZTangent;
LLeftForce = WaterForceAtStaaf(FXMidLeft, FZMidLeft, LRadius);
LRadius = FZMidRight - FZTangent;
LRightForce = WaterForceAtStaaf(FXMidRight, FZMidRight, LRadius);
var result = LLeftForce - LRightForce;
LeftPoreOnSurface = determinePoreOnSurface(FXMidLeft);
RightPoreOnSurface = determinePoreOnSurface(FXMidRight);
return result + (LeftPoreOnSurface - RightPoreOnSurface);
}
private double WaterForceAtStaaf(double AX, double AY, double ARadius)
{
double LYBot;
//double LYBot;
double LPoreBot;
double LForce = 0;
double LPoretop;
LYBot = AY - ARadius;
var LVertical = new TVertical();
LVertical.XCoordinate = AX;
LVertical.ZCoordinate = LYBot;
LVertical.SoilData.AddRange(FSoilData);
LVertical.GammaWater = FInterfaceData.CurrentWaternet.GammaWater;
if (FInterfaceData.CurrentWaternet.PhreaticLineNumber != -1)
{
LVertical.FreaticLine =
FInterfaceData.CurrentWaternet.PLLines[FInterfaceData.CurrentWaternet.PhreaticLineNumber];
}
LVertical.WaterLines.AddRange(FInterfaceData.CurrentWaternet.Waterlines);
LVertical.PLLines = FInterfaceData.CurrentWaternet.PLLines;
LVertical.FillVerticalWithWaterData();
LVertical.SurfaceLine = FInterfaceData.SurfaceLine;
LVertical.InterfaceCalls = FInterfaceData.InterfaceCalls;
LVertical.FillVertical();
var LSoilTop = LVertical.SoilLevel;
var LFreaticLevel = LVertical.Freaticlevel;
LPoreBot = LVertical.DetermineWaternetStress(LYBot);
var LLastZCoor = LYBot;
double LZHead = 0;
double LZ = 0;
double LDz;
while (LVertical.NextWaterLineIsFound(LLastZCoor, ref LZ, ref LZHead))
{
LPoretop = (LZHead - LZ)*FGamWat;
// rectangualr moment if present
var LRectPore = Math.Min(LPoretop, LPoreBot);
LDz = LZ - LLastZCoor;
LForce = LRectPore*LDz;
// triangulair part if present
LForce = 0.5*LDz*Math.Abs(LPoretop - LPoreBot);
LLastZCoor = LZ;
LPoreBot = LPoretop;
}
// Remove moment above surface as this is already teken into acount
if ((LFreaticLevel > LSoilTop))
{
LDz = LFreaticLevel - LSoilTop;
LForce -= 0.5*LDz*LDz*FInterfaceData.CurrentWaternet.GammaWater;
}
return LForce;
}
// Zone calculation
// =======================================================================================================================
// Date ID Modification
// 2010-04-26 Created
// =======================================================================================================================
private void GetRotatedXMiddle()
{
double LVerplaats;
int i;
int j;
double LLength;
bool Ready;
// zoek de afstand waarover de middelpunten verplaatsen
if (FReverseGeom)
{
LVerplaats = Math.Abs(FAngle*FRadRight);
for (i = NSlices - 1; i >= 0; i--)
{
GetSlice(i).IsActive = true;
LLength = 0;
j = i;
do
{
LLength = LLength + GetSlice(j).ArcLength;
Ready = (LLength >= LVerplaats);
if (Ready)
{
// xmidrotated ligt op de plaats van de laatste j
GetSlice(i).XMidRotated = GetSlice(j).xLeft - (LLength - LVerplaats)*(GetSlice(j).xLeft - GetSlice(j).xRight);
}
j -= 1;
} while (!(Ready || (j == -1)));
if (!Ready)
{
GetSlice(i).IsActive = false;
}
}
}
else
{
LVerplaats = Math.Abs(FAngle*FRadLeft);
for (i = 0; i < NSlices; i++)
{
GetSlice(i).IsActive = true;
LLength = 0;
j = i;
do
{
LLength = LLength + GetSlice(j).ArcLength;
Ready = (LLength >= LVerplaats);
if (Ready)
{
// xmidrotated ligt op de plaats van de laatste j
GetSlice(i).XMidRotated = GetSlice(j).xRight - (LLength - LVerplaats)*(GetSlice(j).xRight - GetSlice(j).xLeft);
}
j++;
} while (!(Ready || (j == NSlices)));
if (!Ready)
{
GetSlice(i).IsActive = false;
}
}
}
}
// =======================================================================================================================
// Description: Calculates the DrivingSoilMoment on the left side
//
// Date ID Modification
// 1999-03-29 Best Created
// =======================================================================================================================
private void GetRotatedDrivingMoments(ref double RDriveMoment, ref double LDriveMom)
{
int i;
bool Ready;
LDriveMom = 0;
Ready = false;
i = 0;
while (!Ready && (i < NSlices))
{
Ready = (GetSlice(i).XMidRotated >= FXMidLeft);
if (!Ready)
{
if (GetSlice(i).IsActive)
{
LDriveMom = LDriveMom + GetSlice(i).Weight*(FXMidLeft - GetSlice(i).XMidRotated);
}
}
i++;
}
RDriveMoment = 0;
Ready = false;
i = NSlices - 1;
while (!Ready && (i >= 0))
{
Ready = (GetSlice(i).XMidRotated >= FXMidRight);
if (!Ready)
{
if (GetSlice(i).IsActive)
{
RDriveMoment = RDriveMoment + GetSlice(i).Weight*(FXMidRight - GetSlice(i).XMidRotated);
}
}
i -= 1;
}
}
// =======================================================================================================================
// Description: Calculates the stability factor of a Mlift slide plane
// after rotation.
//
// Date ID Modification
// 1999-03-23 Best Created
// 2006-05-29 Log Check if the driving moment is > 0, otherwise return -1
// =======================================================================================================================
private double GetRotatedStabilityFactor()
{
double result = 0;
int i = 0;
double RDriveMoment = 0;
double LDriveMoment = 0;
double LTotDrivingLeft = 0;
double LTotDrivingRight = 0;
double LDrivingMoment = 0;
double LResisting = 0;
double RResisting = 0;
double Horresisting = 0;
double LResistingMom = 0;
double ArmLeft = 0;
double ArmRight = 0;
// bepaal de nieuwe slice middens
GetRotatedXMiddle();
GetRotatedDrivingMoments(ref RDriveMoment, ref LDriveMoment);
LTotDrivingLeft = LDriveMoment + FWatermomentLeft + FLoadMomentLeft;
LTotDrivingRight = RDriveMoment + FWatermomentRight + FLoadMomentRight;
LResisting = 0;
RResisting = 0;
Horresisting = 0;
// resisting moments
for (i = 0; i < FNSlices; i++)
{
if (GetSlice(i).IsActive)
{
if ((GetSlice(i).XMidRotated < FXMidLeft))
{
LResisting = LResisting + GetSlice(i).ArcLength*GetSlice(i).ShearStress;
}
else if ((GetSlice(i).XMidRotated > FXMidRight))
{
RResisting = RResisting + GetSlice(i).ArcLength*GetSlice(i).ShearStress;
}
else
{
Horresisting = Horresisting + GetSlice(i).ArcLength*GetSlice(i).ShearStress;
}
}
}
// de linker en rechter arm blijven zo
// tgv rotatie zal dat wel iets anders zijn, maar hoe is onbekend
Tslice _wvar1 = FSlices[FLeftNumber];
ArmLeft = FRadLeft - 1/3*(_wvar1.ZTopRight - _wvar1.ZBottomRight);
Tslice _wvar2 = FSlices[FRightNumber];
ArmRight = FRadRight - 1/3*(_wvar2.ZTopLeft - _wvar2.ZBottomLeft);
LResistingMom = Math.Abs(LResisting*FRadLeft)/ArmLeft + Math.Abs(RResisting*FRadRight)/ArmRight + Math.Abs(Horresisting);
LDrivingMoment = LTotDrivingLeft/ArmLeft + LTotDrivingRight/ArmRight;
FTextileMoment = GeotextileMoment(LDrivingMoment);
LResistingMom = LResistingMom + FTextileMoment;
if ((LDrivingMoment > Constants.CGeoAccu))
{
result = LResistingMom/Math.Abs(LDrivingMoment);
}
else
{
result = 1001;
}
return result;
}
// =======================================================================================================================
// Description: ZFunction
// Calculates the
// the variables Ex ed are build as follows :
// in the array Fstrength consists of materials and
// connection to cohesion ore phi.
//
// Date ID Modification
// 1999-09-10 Best Created
// =======================================================================================================================
private void ZFunction(int NTot, double[] X, double theta, ref int IError, ref double Z)
{
double safety = 0;
int iter = 0;
string LStr;
// vul de nieuwe sterkte parameters in en anpassingen
try
{
FStabProbabilisticData.PutAdjustedParametersInMaterial(X);
}
catch (Exception E)
{
LStr = E.Message;
throw E;
}
// Fill soil,stress data of slices (met nieuwe gegevens
// in een apparte routine STSliceFill
try
{
FStabProbabilisticData.AdjustSliceData(FNSlices, ref FSlices);
}
catch (Exception E)
{
LStr = E.Message;
throw E;
}
DetermineDrivingMoments();
// (TotDrivingLeft,TotDrivingRight);
safety = LiftMethode(ref iter);
Z = safety - FStabProbabilisticData.LimitvalueModel;
if ((Z < 900))
{
IError = 0;
}
}
// =======================================================================================================================
// Date ID Modification
// 2008-07-01 Best Created
// =======================================================================================================================
// =======================================================================================================================
// Date ID Modification
// 2008-07-01 Best Created
// =======================================================================================================================
private double GetActiveRadius()
{
double result;
result = Math.Abs(FZActiveCentre - FZTangent);
return result;
}
// =======================================================================================================================
// Date ID Modification
// 2008-07-01 Best Created
// =======================================================================================================================
private double GetPassiveRadius()
{
double result;
result = Math.Abs(FZPassiveCentre - FZTangent);
return result;
}
// =======================================================================================================================
// Date ID Modification
// 2008-07-01 Best Created
// =======================================================================================================================
private double GeotextileMoment(double DrivingMom)
{
double result;
int i;
double Som;
Som = 0;
for (i = FInterfaceData.Geotextiles.GetLowerBound(0); i <= FInterfaceData.Geotextiles.GetUpperBound(0); i++)
{
TGeotextiles _wvar1 = FInterfaceData.Geotextiles[i];
Som = Som + _wvar1.GetGeotextilesMoment(FXActiveCentre, FZActiveCentre, ActiveRadius, DrivingMom);
}
result = Som;
return result;
}
private string TestCenterpointInGeometry(double AXMid, double AZMid)
{
string result;
double LZSurface;
result = "";
LZSurface = MStabDatafunctions.ZBottomAtSurface(AXMid, FInterfaceData.SurfaceLine);
if ((LZSurface >= AZMid))
{
result = "Center point in geometry";
}
return result;
}
// =======================================================================================================================
// Date ID Modification
// 2008-07-01 Best Created
// =======================================================================================================================
private double GetStabilityFactor()
{
double result;
int LIter;
string lstr;
try
{
LIter = 0;
result = Constants.CFMinStart;
PrepareCalculation();
if ((IsSlipPlaneValid()))
{
CreateSlices();
// find slice numbers from start and end of staight line
FindStraightLineSliceNumbers();
// Fill soil,stress data of slices
FillSlices();
// Calculate the driving moments
DetermineDrivingMoments();
// Calculate stability factor
result = LiftMethode(ref LIter);
}
}
catch (Exception E)
{
lstr = E.Message;
throw E;
}
return result;
}
}
// end TUpliftPlane
}