Index: DamEngine/branches/nwo_18.1/pythonScripts/NWODefinition.py =================================================================== diff -u --- DamEngine/branches/nwo_18.1/pythonScripts/NWODefinition.py (revision 0) +++ DamEngine/branches/nwo_18.1/pythonScripts/NWODefinition.py (revision 1980) @@ -0,0 +1,190 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Nov 27 16:14:14 2018 + +@author: noordam +""" +import helpFunctions +import math +import shutil +import DGeoStabParser + +class NWO: + + def __init__(self): + self.localXCoords=[] + self.localYCoords=[] + + self.globalXCoords=[] + self.globalYCoords=[] + + self.NWOXBoundaryType=0 + self.NWOYBoundaryType=0 + self.epsilon=1e-1 + +#============================================================================== +# local coordinates of NWO +#============================================================================== + def CreateCustomNWO(self, h1,n1,h2,n2,b): + x1Loc =0; + y1Loc =0; + + x2Loc = x1Loc+h1 * n1; + y2Loc = y1Loc-h1; + + x3Loc = x2Loc + math.sqrt(b**2-(h2-h1)**2); + y3Loc = y1Loc-h2; + + x4Loc = x3Loc + h2*n2; + y4Loc=y1Loc; + + #bTot = x4Loc-x1Loc; + + self.localXCoords = [x1Loc, x2Loc, x3Loc, x4Loc]; + self.localYCoords = [y1Loc, y2Loc, y3Loc, y4Loc]; + + + def CalculateGlobalCoordinatesCustomNWO(self,newX1Glob,xCoordsOriginalSurfaceLine,yCoordsOriginalSurfaceLine): + x1Glob = newX1Glob + xCoordsLine1, yCoordsLine1, coordIndexLineStart1 = helpFunctions.GetNearestLineAtX(x1Glob,xCoordsOriginalSurfaceLine,yCoordsOriginalSurfaceLine) + y1Glob = helpFunctions.GetYatX(x1Glob,xCoordsLine1,yCoordsLine1) + #============================================================================== + # calculate global coordinates + #============================================================================== + x1Loc=self.localXCoords[0] + x2Loc=self.localXCoords[1] + x3Loc=self.localXCoords[2] + x4Loc=self.localXCoords[3] + + y2Loc=self.localYCoords[1] + y3Loc=self.localYCoords[2] + + bTot = x4Loc-x1Loc; + + + intersectionIsFound = False + ii=0 + while not intersectionIsFound and iix4Glob-self.epsilon: + x3Glob=x4Glob-self.epsilon + + self.globalXCoords = [x1Glob, x2Glob, x3Glob, x4Glob]; + self.globalYCoords = [y1Glob, y2Glob, y3Glob, y4Glob]; + + + def CreateBoundaryAtPleistocene(self,newX1Glob,xCoordsOriginalSurfaceLine,yCoordsOriginalSurfaceLine,fileName): + + x1Glob = newX1Glob + xCoordsLine1, yCoordsLine1, coordIndexLineStart1 = helpFunctions.GetNearestLineAtX(x1Glob,xCoordsOriginalSurfaceLine,yCoordsOriginalSurfaceLine) + y1Glob = helpFunctions.GetYatX(x1Glob,xCoordsLine1,yCoordsLine1) + + # file = r".\R032+20_Bishop_freatische lijn_aangepaste schematisatie.sti" + dgs = DGeoStabParser.DGeoStability(fileName) + + x2Glob= newX1Glob+self.epsilon + y2Glob = dgs.getLayerTop(0, x2Glob) + + x3Glob = xCoordsOriginalSurfaceLine[-1]-self.epsilon + y3Glob = dgs.getLayerTop(0, x3Glob) + + x4Glob = xCoordsOriginalSurfaceLine[-1] + y4Glob = yCoordsOriginalSurfaceLine[-1] + + self.globalXCoords = [x1Glob, x2Glob, x3Glob, x4Glob]; + self.globalYCoords = [y1Glob, y2Glob, y3Glob, y4Glob]; + +# ============================================================================= +# +# +# +# x1Glob = newX1Glob +# xCoordsLine1, yCoordsLine1, coordIndexLineStart1 = helpFunctions.GetNearestLineAtX(x1Glob,xCoordsOriginalSurfaceLine,yCoordsOriginalSurfaceLine) +# y1Glob = helpFunctions.GetYatX(x1Glob,xCoordsLine1,yCoordsLine1) +# #============================================================================== +# # calculate global coordinates +# #============================================================================== +# +# # ============================================================================= +# # Add forbidden line to NWO +# # ============================================================================= +# old_name='D:\\noordam\\Documents\\DAMNWO\\calculations\\test2\\input\\Geometries\\DWP_3.sti' +# new_name= workingDirectory+'Geometries\\DWP_3_'+str(x1Glob)+'.sti' +# shutil.copy(old_name, new_name) +# +# # readAndWriteFiles.AddForbiddenLine(x1Glob,x2Glob,y1Glob,y2Glob,xCoordsNewSurfaceLine,yCoordsNewSurfaceLine,new_name) +# +# #file = r".\R032+20_Bishop_freatische lijn_aangepaste schematisatie.sti" +# dgs = DGeoStabParser.DGeoStability(new_name) +# +# +# +# intersectionIsFound = False +# ii=0 +# while not intersectionIsFound and iix4Glob-epsilon: +# x3Glob=x4Glob-epsilon +# +# globalXCoords = [x1Glob, x2Glob, x3Glob, x4Glob]; +# globalYCoords = [y1Glob, y2Glob, y3Glob, y4Glob]; +# ============================================================================= Index: DamEngine/branches/nwo_18.1/pythonScripts/readAndWriteFiles.py =================================================================== diff -u --- DamEngine/branches/nwo_18.1/pythonScripts/readAndWriteFiles.py (revision 0) +++ DamEngine/branches/nwo_18.1/pythonScripts/readAndWriteFiles.py (revision 1980) @@ -0,0 +1,241 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Feb 20 10:39:51 2018 + +@author: noordam +""" +import csv +import helpFunctions +import re +import sys +import os +#============================================================================== +# read surfaceLine +#============================================================================== +def ReadSurfaceLine(surfaceLineFileName,lineNumber): + with open(surfaceLineFileName,'rb') as csvfile: + reader = csv.reader(csvfile,delimiter=';', quotechar='|') + # print spamreader.row + definitions = next(reader) + for i in range(lineNumber-1): + next(reader) + row = next(reader) + row = filter(None,row) + + originalProfileName =row[0].replace(',','.'); + # geologicalProfile =row.index('geologicalprofile') + +#============================================================================== +# geologicalProfile =row[1].replace(',','.'); +# xGridPoint = row[2].replace(',','.'); +# yGridPoint = row[3].replace(',','.'); +# scenarioClusterID = row[4].replace(',','.'); +#============================================================================== + + additionalInfo = []; + additionalInfo.append(originalProfileName); + +#============================================================================== +# additionalInfo.append(geologicalProfile); +# additionalInfo.append(xGridPoint); +# additionalInfo.append(yGridPoint); +# additionalInfo.append(scenarioClusterID); +#============================================================================== + x1Index = definitions.index('X1') + nColsInputSurfaceLine = len(row) + xCoordsOriginalSurfaceLine = row[x1Index:nColsInputSurfaceLine:3]; + yCoordsOriginalSurfaceLine = row[x1Index+2:nColsInputSurfaceLine:3]; + + xCoordsOriginalSurfaceLine =[float(i.replace(',','.')) for i in xCoordsOriginalSurfaceLine] + yCoordsOriginalSurfaceLine =[float(i.replace(',','.')) for i in yCoordsOriginalSurfaceLine] + + for i in range(1,x1Index): + info = row[i].replace(',','.'); + additionalInfo.append(info); + + + + return xCoordsOriginalSurfaceLine, yCoordsOriginalSurfaceLine, additionalInfo,definitions + +#============================================================================== +# read characteristic points +#============================================================================== +def ReadCharacteristicPoints(charPointsFileName,lineNumber): + with open(charPointsFileName,'rb') as csvfile: + reader = csv.reader(csvfile,delimiter=';', quotechar='|') + # print spamreader.row + definitions = next(reader) + for i in range(lineNumber-1): + next(reader) + row = next(reader) + + nColsInputCharacteristicPoints = len(row) + xCoordsOriginalCharacteristicPoints = row[1:nColsInputCharacteristicPoints-1:3]; + yCoordsOriginalCharacteristicPoints = row[3:nColsInputCharacteristicPoints:3]; + + xCoordsOriginalCharacteristicPoints =[float(i.replace(',','.')) for i in xCoordsOriginalCharacteristicPoints] + yCoordsOriginalCharacteristicPoints =[float(i.replace(',','.')) for i in yCoordsOriginalCharacteristicPoints] + + return xCoordsOriginalCharacteristicPoints, yCoordsOriginalCharacteristicPoints,definitions + + +#============================================================================== +# Read locations, scenarios and segments +#============================================================================== +def ReadOtherCsv(fileName,lineNumber): + with open(fileName,'rb') as csvfile: + reader = csv.reader(csvfile,delimiter=';', quotechar='|') + # print spamreader.row + definitions = next(reader) + for i in range(lineNumber-1): + next(reader) + data = next(reader) + data = [i.replace(',','.') for i in data] + return definitions, data + +#============================================================================== +# write the new surface line to a file +#============================================================================== +def WriteSurfaceLine(nLine, definitions,xCoords,zCoords,additionalInfo,safeFileName): + + newRow =[] + newRow.extend(additionalInfo) + for i in range(len(xCoords)): + newRow.extend([xCoords[i]]) + newRow.extend([0]) + newRow.extend([zCoords[i]]) + + with open(safeFileName,'ab') as csvFile: + writer = csv.writer(csvFile, delimiter =';', quotechar='|') + if (nLine==0): + writer.writerow(definitions) + writer.writerow(newRow) +#============================================================================== +# write characteristicPoints to a new file +#============================================================================== +def WriteCharacteristicPoints(nLine, definitions,xCoords,zCoords,profileName,safeFileName): + #combine all the coordinates in 1 row + newRow =[] + newRow.extend([profileName]) + + for i in range(len(xCoords)): + newRow.extend([xCoords[i]]) + if (xCoords[i]==-1 and zCoords[i]==-1): + newRow.extend([-1]) + else: + newRow.extend([0]) + newRow.extend([zCoords[i]]) + newRow.extend([0]) + + with open(safeFileName,'ab') as csvFile: + writer = csv.writer(csvFile, delimiter =';', quotechar='|') + if (nLine==0): + writer.writerow(definitions) + writer.writerow(newRow) + + return + +# ============================================================================= +# add forbidden line +# ============================================================================= +def AddForbiddenLine(x1,x2,z1,z2,xCoordsSurfaceLine,zCoordsSurfaceLine,fileName): + indexesWithinLine = helpFunctions.FindIndexesInRange(x1,x2,xCoordsSurfaceLine,"surfaceLine") + xCoordsForbiddenLine=[] + zCoordsForbiddenLine=[] + for ii in range(len(indexesWithinLine)): + xCoordsForbiddenLine.append(xCoordsSurfaceLine[indexesWithinLine[ii]]) + zCoordsForbiddenLine.append(zCoordsSurfaceLine[indexesWithinLine[ii]]) + xCoordsForbiddenLine.append(x2) + zCoordsForbiddenLine.append(z2) + + stiFile = open(fileName, 'r') + stiFileLines= stiFile.readlines() + stiFile.close + for ii in range(len(stiFileLines)): + if '[UNIT WEIGHT WATER]' in stiFileLines[ii]: + stiFileLines.insert(ii,'[FORBIDDEN LINES] \n') + ii+=1 + stiFileLines.insert(ii,str(len(xCoordsForbiddenLine)-1) + ' - Forbidden lines \n') + ii+=1 + for jj in range(len(xCoordsForbiddenLine)-1): + stiFileLines.insert(ii,str(jj+1)+ "\n") + ii+=1 + stiFileLines.insert(ii,str(xCoordsForbiddenLine[jj])+' '+str(zCoordsForbiddenLine[jj]) + ' - X1 Y1 \n') + ii+=1 + stiFileLines.insert(ii,str(xCoordsForbiddenLine[jj+1])+' '+str(zCoordsForbiddenLine[jj+1]) + ' - X2 Y2 \n') + ii+=1 + stiFileLines.insert(ii,'[END OF FORBIDDEN LINES] \n') + ii+=1 + break + + #outputStiName = fileName.replace(".sti","_X"+str(x1)+".sti") + outputStiName = fileName + + newStiFile=open(outputStiName,'w') + newStiFile.writelines(stiFileLines) + + test=1+1; +# ============================================================================= +# +# def ReadStiFile(fileName): +# stiFile = open(fileName, 'r') +# stiFileLines= stiFile.readlines() +# stiFile.close +# +# soilLayer=[] +# soilIsAquifer=[] +# aquiferSoils=[] +# topBoundary=[] +# point =[] +# #boundaryTopLayer +# foundSoilInList=False +# foundLayers=False +# +# for ii in range(len(stiFileLines)): +# if (('[SOIL]' in stiFileLines[ii]) or (foundSoilInList == True)): +# if ('[SOIL]' in stiFileLines[ii]): +# foundSoilInList=True +# soilLayer.append(stiFileLines[ii+1]) +# if "SoilIsAquifer" in stiFileLines[ii]: +# soilIsAquifer.append(int(re.search(r'\d+', stiFileLines[ii]).group())) +# if soilIsAquifer[-1] == 1: +# aquiferSoils.append(soilLayer[-1]) +# foundSoilInList = False +# if '[POINTS]' in stiFileLines[ii]: +# nPoints = int(re.search(r'\d+', stiFileLines[ii+1]).group()) +# for jj in range(nPoints): +# point.append([float(x) for x in stiFileLines[ii+jj+2].split()]) +# if '[CURVES]' in stiFileLines[ii]: +# nCurves = int(re.search(r'\d+', stiFileLines[ii+1]).group()) +# for jj in range(nCurves): +# # @@!@!!@!! point.append([float(x) for x in stiFileLines[ii+jj+2].split()]) +# +# if '[LAYERS]' in stiFileLines[ii] or (foundLayers==True): +# foundLayers=True +# for jj in range(len(aquiferSoils)): +# if (aquiferSoils[jj] in stiFileLines[ii+3]): +# topBoundary.append(int(re.search(r'\d+', stiFileLines[ii+3]).group())) +# if (len(aquiferSoils)==len(topBoundary)): +# foundLayers = False +# +# test=1+1 +# +# ============================================================================= + +def WriteOtherCsv(nLine, definitions, data, safeFileName): + with open(safeFileName,'ab') as csvFile: + writer = csv.writer(csvFile, delimiter =';', quotechar='|') + if (nLine==0): + writer.writerow(definitions) + writer.writerow(data) + + +def WritePythonInformationForDamKernel(damExecutableFolder): + with open(damExecutableFolder+'PythonInformation','w') as cSharpFile: + cSharpFile.write(sys.executable+'\n') + cSharpFile.write(os.getcwd()+'\n') + cSharpFile.write('AddForbiddenLine = 1\n') + cSharpFile.write('AddMaxXEntrance = 1\n') + + + \ No newline at end of file Index: DamEngine/branches/nwo_18.1/pythonScripts/GeometryDefinition.py =================================================================== diff -u --- DamEngine/branches/nwo_18.1/pythonScripts/GeometryDefinition.py (revision 0) +++ DamEngine/branches/nwo_18.1/pythonScripts/GeometryDefinition.py (revision 1980) @@ -0,0 +1,156 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Nov 27 16:46:58 2018 + +@author: noordam +""" +import readAndWriteFiles +import helpFunctions + +class SurfaceLine: + + def __init__(self): + self.xCoordsOriginalSurfaceLine=[] + self.yCoordsOriginalSurfaceLine=[] + self.xCoordsNewSurfaceLine=[] + self.yCoordsNewSurfaceLine=[] + + self.info=[] + self.surfaceLineDefinitions =[] + + def ReadSurfaceLine(self,surfaceLineFileName,lineInOriginalFile) : + self.xCoordsOriginalSurfaceLine, self.yCoordsOriginalSurfaceLine, self.info, self.surfaceLineDefinitions = readAndWriteFiles.ReadSurfaceLine(surfaceLineFileName,lineInOriginalFile) + + + def WriteSurfaceLine(self,jj,safeFileNameNewSurfaceLine): + readAndWriteFiles.WriteSurfaceLine(jj,self.surfaceLineDefinitions,self.xCoordsNewSurfaceLine,self.yCoordsNewSurfaceLine,self.info,safeFileNameNewSurfaceLine) + #readAndWriteFiles.WriteCharacteristicPoints(jj,list(reversed(newCharaterisitcDefinitions)),list(reversed(xCoordsNewCharacteristicPoints)),list(reversed(yCoordsNewCharacteristicPoints)),profileName,safeFileNameCharacteristicPoints) + + + def CreateNewSurfaceLine(self,nwo): + + x1Glob = nwo.globalXCoords[0] + x2Glob = nwo.globalXCoords[1] + x3Glob = nwo.globalXCoords[2] + x4Glob = nwo.globalXCoords[3] + + y1Glob = nwo.globalYCoords[0] + y2Glob = nwo.globalYCoords[1] + y3Glob = nwo.globalYCoords[2] + y4Glob = nwo.globalYCoords[3] + + #============================================================================== + # Create new surfaceLine + #============================================================================== + self.xCoordsNewSurfaceLine = self.xCoordsOriginalSurfaceLine[:] + self.yCoordsNewSurfaceLine = self.yCoordsOriginalSurfaceLine[:] + + # replace the y coordinates of the existing points + indexesWithinNWO = helpFunctions.FindIndexesInRange(x1Glob,x4Glob,self.xCoordsOriginalSurfaceLine,"surfaceLine") + for ii in range(len(indexesWithinNWO)): + if self.xCoordsNewSurfaceLine[indexesWithinNWO[ii]]=0: + if "Fmin" not in textInFile[j]: + j=j-1; + continue + else: + fsIsFound = True + [float(s) for s in textInFile[j].split() if s.isdigit()] + fMin =float(s); + allfMin.append(fMin) + + \ No newline at end of file Index: DamEngine/branches/nwo_18.1/pythonScripts/DGeoStabParser.py =================================================================== diff -u --- DamEngine/branches/nwo_18.1/pythonScripts/DGeoStabParser.py (revision 0) +++ DamEngine/branches/nwo_18.1/pythonScripts/DGeoStabParser.py (revision 1980) @@ -0,0 +1,235 @@ +import matplotlib.pyplot as plt +import numpy as np + +class Soil: + + def __init__(self): + self.SoilName = "" + self.SoilColor = 00000000 + self.SoilSoilType = 0 + self.SoilUseSoilType = 0 + self.SoilExcessPorePressure = 0.00 + self.SoilPorePressureFactor = 0.00 + self.SoilGamDry = 0.00 + self.SoilGamWet = 0.00 + self.SoilRestSlope = 0 + self.SoilCohesion = 0.00 + self.SoilPhi = 0.00 + self.SoilDilatancy = 0.00 + self.SoilCuTop = 0.00 + self.SoilCuBottom = 0.00 + self.SoilCuGradient = 0.00 + self.SoilStressTableName = "" + self.SoilBondStressTableName = "" + self.SoilMatStrengthType = 0 + self.SoilProbInputValues = 0 + self.SoilRatioCuPc = 0.0 + self.SoilStrengthIncreaseExponent = 0.0 + self.SoilPOPTop = 0.00 + self.SoilPOPBottom = 0.00 + self.SoilAdditionalFactorLEM = 0.00 + self.SoilPc = 0.00E+00 + self.SoilPOP = 0.00 + self.SoilRheologicalCoefficient = 0.00 + self.SoilXCoorSoilPc = 0.000 + self.SoilYCoorSoilPc = 0.000 + self.SoilIsPopCalculated = 0 + self.SoilIsOCRCalculated = 0 + self.SoilIsAquifer = 0 + self.SoilUseProbDefaults = 0 + self.SoilStdCohesion = 0.00 + self.SoilStdPhi = 0.00 + self.SoilStdRatioCuPc = 0.00 + self.SoilStdRatioCuPcPassive = 0.00 + self.SoilStdRatioCuPcActive = 0.00 + self.SoilStdCu = 0.00 + self.SoilStdCuTop = 0.00 + self.SoilStdCuGradient = 0.00 + self.SoilStdPn = 0.20 + self.SoilDistCohesion = 0 + self.SoilDistPhi = 0 + self.SoilDistStressTable = 0 + self.SoilDistRatioCuPc = 0 + self.SoilDistRatioCuPcPassive = 0 + self.SoilDistRatioCuPcActive = 0 + self.SoilDistCu = 0 + self.SoilDistCuTop = 0 + self.SoilDistCuGradient = 0 + self.SoilDistPn = 0 + self.SoilCorrelationCPhi = 0.00 + self.SoilRatioCuPcPassive = 0.00 + self.SoilRatioCuPcActive = 0.00 + self.SoilCuPassiveTop = 0.00 + self.SoilCuPassiveBottom = 0.00 + self.SoilCuActiveTop = 0.00 + self.SoilCuActiveBottom = 0.00 + self.SoilUniformRatioCuPc = 0 + self.SoilUniformCu = 0 + self.SoilDesignPartialCohesion = 0.00 + self.SoilDesignStdCohesion = 0.00 + self.SoilDesignPartialPhi = 0.00 + self.SoilDesignStdPhi = 0.00 + self.SoilDesignPartialStressTable = 0.00 + self.SoilDesignStdStressTable = 0.00 + self.SoilDesignPartialRatioCuPc = 0.00 + self.SoilDesignStdRatioCuPc = 0.00 + self.SoilDesignPartialCu = 0.00 + self.SoilDesignStdCu = 0.00 + self.SoilDesignPartialPOP = 0.00 + self.SoilDesignStdPOP = 0.00 + self.SoilDesignPartialRRatio = 0.00 + self.SoilDesignStdRRatio = 0.00 + self.SoilSoilGroup = 0 + self.SoilStdPOP = 0.00 + self.SoilDistPOP = 0 + self.SoilHorFluctScaleCoh = 0.00 + self.SoilVertFluctScaleCoh = 0.00 + self.SoilNumberOfTestsCoh = 0 + self.SoilVarianceRatioCoh = 0.00 + self.SoilHorFluctScalePhi = 0.00 + self.SoilVertFluctScalePhi = 0.00 + self.SoilNumberOfTestsPhi = 0 + self.SoilVarianceRatioPhi = 0.00 + self.SoilRRatio = 0.0000000 + self.SoilDistCu = 0 + self.SoilDistCuTop = 0 + self.SoilDistCuGradient = 0 + ## self.MPMMat = Materials.Materials() + + def read_soil(self, file): + line = file.readline() + self.SoilName = line.strip() + line = file.readline() + while line.strip() != "[END OF SOIL]": + cells = line.strip().split('=') + if len(cells) == 2: + if cells[1] != '': + exec("self." + line) + elif len(cells) > 2: + exec("self." + cells[0] + "=" + "\"" + " ".join(cells[1:]) + "\"") + line = file.readline() + +class Geometry: + + def __init__(self): + self.Points = [[-1, -1, -1]] + self.Curves = [[[], 0]] + self.Bounds = [] + self.Layers = [] + + + def read_points(self, file): + file.readline() + line = file.readline() + ptName = 1 + while line.strip() != "[END OF POINTS]": + cells = line.split() + self.Points.append([float(cells[1]), float(cells[2]), ptName]) + ptName += 1 + line = file.readline() + + def read_curves(self, file): + file.readline() + + file.readline() + line = file.readline() + + lineName = 1 + while line.strip() != "[END OF CURVES]": + line = file.readline() + cells = line.split() + self.Curves.append([[int(val) for val in cells], lineName]) + lineName += 1 + line = file.readline() + file.readline() + + def read_boundaries(self, file): + # number of boundaries + noBounds = int(file.readline().split()[0]) + for bound in range(noBounds): + file.readline() + noVals = int(file.readline().split()[0]) + noLines = int(noVals / 10) + if noVals % 10 != 0: + noLines += 1 + vals = [] + for i in range(noLines): + line = file.readline().strip() + cells = line.split() + for cell in cells: + vals.append(int(cell)) + self.Bounds.append(vals) + + + def read_layers(self, file): + # number of layers + noLayers = int(file.readline().split()[0]) + + for _ in range(noLayers): + file.readline() + soil = file.readline().strip() + piezoTop = int(file.readline().split()[0]) + piezoBottom = int(file.readline().split()[0]) + topBound = int(file.readline().split()[0]) + bottomBound = int(file.readline().split()[0]) + self.Layers.append([topBound, bottomBound, soil, piezoTop, piezoBottom]) + +class DGeoStability: + + def __init__(self, path): + self.file = open(path, 'r') + self.Geometry = Geometry() + self.readfile() + + def readfile(self): + line = self.file.readline() + while line.strip() != "[END OF INPUT FILE]": + + if line.strip() == "[POINTS]": + self.Geometry.read_points(self.file) + + elif line.strip() == "[CURVES]": + self.Geometry.read_curves(self.file) + + elif line.strip() == "[BOUNDARIES]": + self.Geometry.read_boundaries(self.file) + + elif line.strip() == "[LAYERS]": + self.Geometry.read_layers(self.file) + + line = self.file.readline() + + def getLayerTop(self, layer, x): + topLine = self.Geometry.Layers[layer][0] + curves = self.Geometry.Bounds[topLine] + points = [self.Geometry.Curves[curves[0]][0][0]] + for curve in curves: + points.append(self.Geometry.Curves[curve][0][1]) + + x_pts = [] + z_pts = [] + for point in points: + x_pts.append(self.Geometry.Points[point][0]) + z_pts.append(self.Geometry.Points[point][1]) + + minX = min(x_pts) + maxX = max(x_pts) + + if x > maxX: + #raise Exception("x-value - outside bounds (Right Boundary =" + str(maxX) + ")") + print ("x-value - outside bounds (Left Boundary =" + str(minX) + ")") + x=maxX-0.1 + elif x < minX: + # raise Exception("x-value - outside bounds (Left Boundary =" + str(minX) + ")") + print ("x-value - outside bounds (Left Boundary =" + str(minX) + ")") + x=minX+0.1 + + return np.interp(x, x_pts, z_pts) + +if __name__ == "__main__": + + file = r".\R032+20_Bishop_freatische lijn_aangepaste schematisatie.sti" + dgs = DGeoStability(file) + z = dgs.getLayerTop(4, 155.0) # DGeoStability( layerNo, x_val) + # print (z) + Index: DamEngine/branches/nwo_18.1/pythonScripts/ontgrondingsProfiel.py =================================================================== diff -u --- DamEngine/branches/nwo_18.1/pythonScripts/ontgrondingsProfiel.py (revision 0) +++ DamEngine/branches/nwo_18.1/pythonScripts/ontgrondingsProfiel.py (revision 1980) @@ -0,0 +1,198 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Feb 19 15:44:34 2018 + +@author: noordam +""" + +import matplotlib.pyplot as plot +import readAndWriteFiles +import os +import NWODefinition +import GeometryDefinition +import sys + + + +epsilon=1e-1 # difference in x coordinate between bottom and top point of NWO, when epsilon is too small, problems in d-geo stability occur. +#============================================================================== +# input +#============================================================================== + +# define the dimensions of the NWO +#============================================================================== +# || 1 4 || +# h1 || \ / || h2 +# || \1:n1 1:n2/ || +# || \ / || +# 2----------3 +# b +#============================================================================== +# dimensions of the NWO + + + +h1 = 4.5; +n1 =0; +h2 = 4.5; +n2 = 5; +b=5; + + +xIni =58.0; # initial location of the NWO point 1 +dX=2; # x-distance of NWO point 1 between calculations +nSteps =10; # number of calculations + +boundaryType=4 +# define the profile from the original csv files which has to be copied +lineInOriginalFile=3; + +workingDirectory ='D:\\noordam\\Documents\\DAMNWO\\calculations\\test2\\input\\' +damExecutableFolder = 'D:\\noordam\\Documents\\nwoUI_18.1\\lib\\DamEngine\\' + +#============================================================================== +# define the used paths +#============================================================================== +if not os.path.isdir(workingDirectory+'originalCSV\\'): + os.makedirs(workingDirectory+'originalCSV\\') +if not os.path.isdir(workingDirectory+'newCSV\\'): + os.makedirs(workingDirectory+'newCSV\\') + +originalFileLocation = workingDirectory +'originalCSV\\' +newFileLocation = workingDirectory +'newCSV\\' +surfaceLineFileName = originalFileLocation+'surfacelines.csv' +charPointsFileName = originalFileLocation+'characteristicpoints.csv' + +locationsFileName = originalFileLocation+'locations.csv' +scenariosFileName = originalFileLocation+'scenarios.csv' +segmentsFileName = originalFileLocation+'segments.csv' + +safeFileNameCharacteristicPoints = newFileLocation+'characteristicpoints.csv' +safeFileNameNewSurfaceLine = newFileLocation+'surfacelines.csv' + +safeFileNameLocations = newFileLocation+ 'locations.csv' +safeFileNameScenarios = newFileLocation+ 'scenarios.csv' +safeFileNameSegments = newFileLocation+ 'segments.csv' + +#============================================================================== +# remove the to be created file if it already exists +#============================================================================== + +if os.path.isfile(safeFileNameCharacteristicPoints): + os.remove(safeFileNameCharacteristicPoints) +if os.path.isfile(safeFileNameNewSurfaceLine): + os.remove(safeFileNameNewSurfaceLine) + +if os.path.isfile(safeFileNameLocations): + os.remove(safeFileNameLocations) +if os.path.isfile(safeFileNameScenarios): + os.remove(safeFileNameScenarios) +if os.path.isfile(safeFileNameSegments): + os.remove(safeFileNameSegments) + + +#============================================================================== +# local coordinates of NWO +#============================================================================== +nwo = NWODefinition.NWO() +nwo.CreateCustomNWO(h1,n1,h2,n2,b) + +#============================================================================== +# readSurfaceLine and characteristic points +#============================================================================== + +surfaceLine=GeometryDefinition.SurfaceLine() +characteristicPoints=GeometryDefinition.CharacteristicPoints() + +surfaceLine.ReadSurfaceLine(surfaceLineFileName,lineInOriginalFile) +characteristicPoints.ReadCharacteristicPoints(charPointsFileName,lineInOriginalFile) + + +originalProfileName = surfaceLine.info[0] +newX1Glob = xIni + +if (boundaryType ==4): + ii= characteristicPoints.characteristicDefinitions.index('X_Teen dijk binnenwaarts') + newX1Glob = characteristicPoints.xCoordsOriginalCharacteristicPoints[(ii-1)/3] +#============================================================================== +# read the other csv files, locations scenarios and segments +#============================================================================== +locationsDefinitions, locations = readAndWriteFiles.ReadOtherCsv(locationsFileName,lineInOriginalFile) +scenariosDefinitions, scenarios = readAndWriteFiles.ReadOtherCsv(scenariosFileName,lineInOriginalFile) +segmentsDefinitions, segments = readAndWriteFiles.ReadOtherCsv(segmentsFileName,lineInOriginalFile) + +#============================================================================== +# loop over the amount of steps +#============================================================================== + +for jj in range(nSteps): + #============================================================================== + # Find Global coordinate of point 1 + #============================================================================== + + + #nwo.CalculateGlobalCoordinatesCustomNWO(newX1Glob,surfaceLine.xCoordsOriginalSurfaceLine,surfaceLine.yCoordsOriginalSurfaceLine) + fileName = workingDirectory+'Geometries\\DWP_3.sti' + nwo.CreateBoundaryAtPleistocene(newX1Glob,surfaceLine.xCoordsOriginalSurfaceLine,surfaceLine.yCoordsOriginalSurfaceLine,fileName) + + surfaceLine.CreateNewSurfaceLine(nwo) + + ### todo: work on max boundary nwo + characteristicPoints.CreateNewCharacteristicPoints(nwo) + + + + #readAndWriteFiles.ReadStiFile(new_name) + + #============================================================================== + # write new surface line and characteristic points + #============================================================================== + profileName = originalProfileName + profileName = profileName+'xDitchDike'+ '%.2f' % nwo.globalXCoords[0] + profileName = profileName.replace(".", "_") + + surfaceLine.info[0] = profileName + + surfaceLine.WriteSurfaceLine(jj,safeFileNameNewSurfaceLine) + characteristicPoints.WriteCharacteristicPoints(jj,profileName,safeFileNameCharacteristicPoints) + + + #============================================================================== + # write the new locations, scenarios and segments + #============================================================================== + locations[0] = profileName + locations[1] = profileName + locations[2] = profileName + + scenarios[0] = profileName + + segments[0] = profileName + + readAndWriteFiles.WriteOtherCsv(jj, locationsDefinitions, locations, safeFileNameLocations) + readAndWriteFiles.WriteOtherCsv(jj, scenariosDefinitions, scenarios, safeFileNameScenarios) + readAndWriteFiles.WriteOtherCsv(jj, segmentsDefinitions, segments, safeFileNameSegments) + + + #============================================================================== + # update position + #============================================================================== + newX1Glob=nwo.globalXCoords[0]+dX; +#============================================================================== +# plots +#============================================================================== +plot.figure(1) +plot.plot(nwo.localXCoords, nwo.localYCoords) + +plot.figure(2) +plot.plot(nwo.globalXCoords, nwo.globalYCoords) +plot.plot(surfaceLine.xCoordsOriginalSurfaceLine, surfaceLine.yCoordsOriginalSurfaceLine) + +plot.figure(3) +plot.plot(surfaceLine.xCoordsNewSurfaceLine, surfaceLine.yCoordsNewSurfaceLine) +plot.plot(characteristicPoints.xCoordsNewCharacteristicPoints, characteristicPoints.yCoordsNewCharacteristicPoints,'o') + + +readAndWriteFiles.WritePythonInformationForDamKernel(damExecutableFolder) + + +print( "calculation succeeded") \ No newline at end of file Index: DamEngine/branches/nwo_18.1/pythonScripts/helpFunctions.py =================================================================== diff -u --- DamEngine/branches/nwo_18.1/pythonScripts/helpFunctions.py (revision 0) +++ DamEngine/branches/nwo_18.1/pythonScripts/helpFunctions.py (revision 1980) @@ -0,0 +1,121 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Feb 19 17:59:31 2018 + +@author: noordam +""" + +import math + +epsilon =1e-17 +#============================================================================== +# compute the intersection between a line and a circle +#============================================================================== +def computeCircleLineIntersection(xA,yA,xB, yB, xC, yC, r): + + # compute the euclidean distance between A and B + lAB =math.sqrt((xB-xA)**2.0+(yB-yA)**2.0); + + # compute the direction vector D from A to B + dX = (xB - xA) /lAB; + dY = (yB - yA) /lAB; + + # Now the line equation is x = Dx*t + Ax, y = Dy*t + Ay with 0 <= t <= 1. + # compute the value t of the closest point to the circle center (Cx, Cy) + t = dX*(xC -xA) + dY*(yC-yA); + + # This is the projection of C on the line from A to B. + # compute the coordinates of the point E on line and closest to C + xE = t*dX+xA; + yE = t*dY + yA; + + # compute the euclidean distance from E to C + lEC = math.sqrt((xE-xC)**2.0+(yE-yC)**2.0); + + # test if the line intersects the circle + if lEC < r: + # compute distance from t to circle intersection point + dt = math.sqrt( r**2.0 - lEC**2.0) + + # compute first intersection point + # xF = (t-dt)*dX + xA + # yF = (t-dt)*dY + yA + + # compute second intersection point + xG = (t+dt)*dX + xA + yG = (t+dt)*dY + yA + # print xF + # print yF + if xG> xA and xG<=xB: + return(xG,yG) + else: + return False + + # else test if the line is tangent to circle + elif( lEC == r ): + # tangent point to circle is E + print ("only 1 intersection") + return False + else: + # line doesn't touch circle + print ("no intersections") + return False +#============================================================================== +# get the part of the surface line on which x is located +#============================================================================== +def GetNearestLineAtX(x,xCoordsSurfaceLine,yCoordsSurfaceLine): + ii =0; + if x>=xCoordsSurfaceLine[len(xCoordsSurfaceLine)-1]: + return False + + xLinePoint1=0.0; + while x> xCoordsSurfaceLine[ii] and iixCoordsSurfaceLine[ii] and ii=xStart-epsilon: + indexesInRange.extend([ii]) + ii=ii+1; + else: + while ii=xStart-epsilon and xEnd>xCoordsSurfaceLine[ii]: + indexesInRange.extend([ii]) + ii=ii+1; + + + return indexesInRange + + \ No newline at end of file Index: DamEngine/branches/nwo_18.1/pythonScripts/computeCircleLineIntersection.py =================================================================== diff -u --- DamEngine/branches/nwo_18.1/pythonScripts/computeCircleLineIntersection.py (revision 0) +++ DamEngine/branches/nwo_18.1/pythonScripts/computeCircleLineIntersection.py (revision 1980) @@ -0,0 +1,65 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Feb 19 17:59:31 2018 + +@author: noordam +""" + +import math + +def computeCircleLineIntersection(xA,yA,xB, yB, xC, yC, r): +#============================================================================== +# xA = 2; # line point 1 +# yA = 2; +# +# yB = 3; # line point 2 +# xB = 5; +# +# xC = 1; # circle coordinates +# yC = 2; +# r = 2; # radius +#============================================================================== + + # compute the euclidean distance between A and B + lAB =math.sqrt((xB-xA)**2+(yB-yA)**2); + + # compute the direction vector D from A to B + dX = (xB - xA) /lAB; + dY = (yB - yA) /lAB; + + # Now the line equation is x = Dx*t + Ax, y = Dy*t + Ay with 0 <= t <= 1. + # compute the value t of the closest point to the circle center (Cx, Cy) + t = dX*(xC -xA) + dY*(yC-yA); + + # This is the projection of C on the line from A to B. + # compute the coordinates of the point E on line and closest to C + xE = t*dX+xA; + yE = t*dY + yA; + + # compute the euclidean distance from E to C + lEC = math.sqrt((xE-xC)**2+(yE-yC)**2); + + # test if the line intersects the circle + if lEC < r: + # compute distance from t to circle intersection point + dt = math.sqrt( r**2 - lEC**2) + + # compute first intersection point + # xF = (t-dt)*dX + xA + # yF = (t-dt)*dY + yA + + # compute second intersection point + xG = (t+dt)*dX + xA + yG = (t+dt)*dY + yA + # print xF + # print yF + + return(xG,yG) + + # else test if the line is tangent to circle + elif( lEC == r ): + # tangent point to circle is E + return False + else: + # line doesn't touch circle + return False \ No newline at end of file Index: DamEngine/branches/nwo_18.1/pythonScripts/NWO.py =================================================================== diff -u --- DamEngine/branches/nwo_18.1/pythonScripts/NWO.py (revision 0) +++ DamEngine/branches/nwo_18.1/pythonScripts/NWO.py (revision 1980) @@ -0,0 +1,21 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Nov 27 16:14:14 2018 + +@author: noordam +""" + +class NWO: + + def __init__(self): + self.localXCoords=[] + self.localYCoords=[] + + self.globalXCoords=[] + self.globalYCoords=[] + + self.NWOXBoundaryType=0 + self.NWOYBoundaryType=0 + + + \ No newline at end of file Index: DamEngine/branches/nwo_18.1/pythonScripts/DefineXEntryMax.py =================================================================== diff -u --- DamEngine/branches/nwo_18.1/pythonScripts/DefineXEntryMax.py (revision 0) +++ DamEngine/branches/nwo_18.1/pythonScripts/DefineXEntryMax.py (revision 1980) @@ -0,0 +1,36 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Dec 5 09:27:07 2018 + +@author: noordam +""" +import sys + +def DefineXEntryMax(xCrest,fileName): + stiFile = open(fileName, 'r') + stiFileLines= stiFile.readlines() + stiFile.close + ii=0 + while "[END OF SLIP CIRCLE SELECTION]" not in stiFileLines[ii]: + if 'IsMaxXEntryUsed=0' in stiFileLines[ii]: + stiFileLines[ii] = 'IsMaxXEntryUsed=1 \n' + if 'XEntryMax' in stiFileLines[ii]: + stiFileLines[ii] = 'XEntryMax='+str(xCrest)+' \n' + ii=ii+1 + + + outputStiName = fileName + + newStiFile=open(outputStiName,'w') + newStiFile.writelines(stiFileLines) + newStiFile.close(); + + +xCrest=sys.argv[1] +fileName = sys.argv[2]; + + +#xCrest = str(41.9) +#fileName ="'D:\noordam\Documents\DAMNWO\calculations\test2\input\DAMTestNewInput.Calc\Stability\Bishop\Loc(DWP_3xDitchDike74_69)_Sce(1)_Pro(DWP_3).sti' + +DefineXEntryMax(xCrest,fileName) \ No newline at end of file