From 6171efea6b47507f7f8c21ade7b92306d789c4a9 Mon Sep 17 00:00:00 2001 From: cslotboom <50510141+cslotboom@users.noreply.github.com> Date: Sat, 10 Feb 2024 15:19:30 -0800 Subject: [PATCH] bump version to 1.3.0 --- src/Changelog.txt | 11 +- src/analysis.py | 441 ------------------------------------ src/builder.py | 7 +- src/postprocess/__init__.py | 8 +- src/postprocess/parse.py | 98 ++++++++ src/postprocess/plot.py | 241 +++++++++----------- src/postprocess/poi.py | 297 +++++++++++++++++++++++- src/requirements.txt | 2 + 8 files changed, 516 insertions(+), 589 deletions(-) delete mode 100644 src/analysis.py create mode 100644 src/postprocess/parse.py diff --git a/src/Changelog.txt b/src/Changelog.txt index 7e503f0..13bc4ca 100644 --- a/src/Changelog.txt +++ b/src/Changelog.txt @@ -1,11 +1,10 @@ -1.1 release +1.3 release =========== -Loads: -Linear loads are needed for snow drifts and retaining wall design! -They're been added in the latest patch, along with t -- Added Linear Loads! -- Added (vertical) linear loads to the beam diagram. +Diagram labels: +- Added a dependency for the libary textalloc +- Added a label feature to plots, which can be used to find points of interest. + Fixities: Working with a list of fixities is a little awkward. diff --git a/src/analysis.py b/src/analysis.py deleted file mode 100644 index 5d79d7c..0000000 --- a/src/analysis.py +++ /dev/null @@ -1,441 +0,0 @@ -import openseespy.opensees as op -import numpy as np -import matplotlib.pyplot as plt - -import planesections.builder as bb - -from abc import ABC, abstractmethod - - -class OutputRecorder(ABC): - Nnodes:int - nodeID0:float - nodeIDEnd:int - ndf:int - node:list - - @abstractmethod - def getEleInteralForce(): - pass - - -class OutputRecorderOpenSees(OutputRecorder): - """ - An interface that can be used to get beam internal forces for each node - in the model. - When called on a beam, it will get all internal forces for that beam. - Information at each node in the model is stored in the beam. - The recorder is only is not instantiated at the time of recording. - - Parameters - ---------- - beam : planesections Beam2D - The beam whose data is being recorded. - - """ - def __init__(self, beam:bb.Beam): - - - self.Nnodes = beam.Nnodes - self.nodeID0 = 1 - self.nodeIDEnd = self.Nnodes - self.ndf = beam._ndf - - for ii, node in enumerate(beam.nodes): - ID = node.ID - node.disps = np.array(op.nodeDisp(ID)) - node.rFrc = np.array(op.nodeReaction(ID)) - node.Fint = self.getEleInteralForce(ID) - - def getEleInteralForce(self, nodID:int): - """ - Gets the internal force at the left and right side of a node. - The left and right side forces represent internal force at either side - of a section cut. - - """ - ndf = self.ndf - Fint = np.zeros(ndf*2) - if nodID == self.nodeID0: # Left most node - eleR = nodID - # 0 is used to so that the plot "closes", i.e. starts at zero the goes up - Fint[:ndf] = 0 # Left side forces - Fint[ndf:] = op.eleForce(eleR)[:ndf] # Right side forces - - #Direct Check, this is scary. - elif nodID == self.nodeIDEnd: # right side node - eleL = nodID - 1 - Fint[:ndf] = -np.array(op.eleForce(eleL)[ndf:]) # Right side forces - Fint[ndf:] = 0 # Left side forces - else: # center nodes - eleL = nodID - 1 - eleR = nodID - - Fint[:ndf] = -np.array(op.eleForce(eleL)[ndf:]) # left entries - Fint[ndf:] = np.array(op.eleForce(eleR)[:ndf]) # right entries - - return Fint - - -class OutputRecorder2D(OutputRecorderOpenSees): - - def __post_init__(self): - print('OutputRecorder2D is depcricated and will be removed in a future release. Use OutputRecorder instead') - - -class BeamAnalyzer: - pass - - -class OpenSeesAnalyzer2D(BeamAnalyzer): - """ - This class is used to can be used to create and run an analysis of an - input 2D beam using OpenSeesPy. The nodes, elements, sections, and - forces for the beam are defined in the analysis model - - Note, nodes and elements will both start at 0 instead of 1. - - Parameters - ---------- - beam : planesections Beam2D - The beam whose data is being recorded. - recorder : planesections Recorder - The recorder to use for the output beam. - geomTransform: str, optional - The OpenSees Geometry transform to use. Can be "Linear" or "PDelta" - clearOld : bool, optional - A flag that can be used to turn on or off clearing the old analysis - when the beam is created. - There are some very niche cases where users may want to have mutiple - beams at once in the OpenSees model. - However, this should remain true for nearly all analyses. - Do not turn on unless you know what you're doing. - - """ - def __init__(self, beam2D:bb.Beam, recorder = OutputRecorderOpenSees, - geomTransform = 'Linear', clearOld = True): - - self.beam:bb.Beam = beam2D - self._checkBeam(beam2D) - - self.Nnode = beam2D.Nnodes - self.Nele = self.Nnode - 1 - self.recorder = recorder - self.geomTransform = geomTransform - self.clearOld = clearOld - - def _checkBeam(self, beam2D): - if not beam2D._dimension: - raise Exception("The beam has no dimension, something terible has happened.") - if beam2D._dimension != '2D': - raise Exception(f"The beam has dimension of type {beam2D._dimension}, it should have type '2D'") - - def initModel(self, clearOld = True): - """ - Initializes the model. - - Parameters - ---------- - clearOld : bool, optional - A flag that can be used to turn on or off clearing the old analysis - when the beam is created. - There are some very niche cases where users may want to have mutiple - beams at once in the OpenSees model. - However, this should remain true for nearly all analyses. - Do not turn on unless you know what you're doing. - """ - - if clearOld: - op.wipe() - op.model('Basic' , '-ndm', 2) - op.geomTransf(self.geomTransform, 1) - - - def buildNodes(self): - """ - Adds each node in the beam to the OpenSeesPy model, and assigns - that node a fixity. - """ - - for node in self.beam.nodes: - op.node(int(node.ID), float(node.x), 0.) - - # OpenSees is very finicky with these inputs, int them for saftey. - f1, f2, f3 = node.fixity.fixityValues - op.fix(node.ID, int(f1), int(f2), int(f3)) - - def buildEulerBeams(self): - """ - Creates an elastic Euler beam between each node in the model. - """ - beam = self.beam - matPropreties = beam.getMaterialPropreties() - for ii in range(self.Nele): - ID = ii + 1 - eleID = int(ID) - Ni = int(ID) - Nj = int(ID + 1) - op.element(beam.EleType, eleID , Ni, Nj , *matPropreties, 1) - - - def _buildPointLoads(self, pointLoads): - for load in pointLoads: - op.load(int(load.nodeID), *load.P) - - - def buildPointLoads(self): - """ - Applies point loads to the appropriate nodes in the model. - """ - op.timeSeries('Linear',1) - op.pattern('Plain', 1, 1) - self._buildPointLoads(self.beam.pointLoads) - - def buildAnalysisPropreties(self): - """ - Typical openSeesPy propreties that should work for any linear beam. - A linear algorithm is used because there is no nonlienarity in the beam. - """ - # op.constraints("Transformation") - op.constraints("Lagrange") - op.numberer("Plain") - op.system('BandGeneral') - op.test('NormDispIncr', 1.*10**-8, 40, 0 , 2) - # op.algorithm('Newton') - op.algorithm('Linear') - op.integrator('LoadControl',1.) - op.analysis('Static') - - def analyze(self): - """ - Analyzes the model once and records outputs. - """ - ok = op.analyze(1) - op.reactions() - return ok - - def buildEleLoads(self): - """ - Applies element loads to the appropriate elements in the model. - """ - op.timeSeries('Linear',2) - op.pattern('Plain', 2, 2) - - for eleload in self.beam.eleLoads: - N1 = self.beam._findNode(eleload.x1) + 1 - N2 = self.beam._findNode(eleload.x2) + 1 - build = self._selectLoad(eleload) - build([N1, N2], eleload) - - def _selectLoad(self, eleload): - if isinstance(eleload, bb.EleLoadDist): - return self._buildDistLoad - if isinstance(eleload, bb.EleLoadLinear): - return self._buildLinLoad - - def _buildDistLoad(self, Nodes:list[int], eleload:bb.EleLoadLinear): - load = eleload.P - N1, N2 = Nodes[0], Nodes[1] - for ii in range(N1, N2): - op.eleLoad('-ele', int(ii), - '-type', '-beamUniform', load[1], load[0]) - - def _buildLinLoad(self, Nodes:list[int], eleload:bb.EleLoadLinear): - load = eleload.P - N1, N2 = Nodes[0], Nodes[1] - for ii in range(N1, N2): - Node1 = self.beam.nodes[ii - 1] - Node2 = self.beam.nodes[ii] - qx1, qx2 = eleload.getLoadComponents(Node1.x, Node2.x, load[0]) - qy1, qy2 = eleload.getLoadComponents(Node1.x, Node2.x, load[1]) - - aOverL = 0. - bOverL = 1. - op.eleLoad('-ele',int(ii), - '-type','beamUniform', - qy1, qx1, aOverL, bOverL, qy2, qx2) - - def runAnalysis(self, recordOutput = True): - """ - Makes and analyzes the beam in OpenSees. - - Returns - ------- - None. - - """ - - self.initModel(self.clearOld) - self.buildNodes() - self.buildEulerBeams() - self.buildPointLoads() - self.buildEleLoads() - self.buildAnalysisPropreties() - self.analyze() - - if recordOutput == True: - self.recorder(self.beam) - - -class OpenSeesAnalyzer3D(BeamAnalyzer): - """ - This class is used to can be used to create and run an analysis of an - input 2D beam using OpenSeesPy. The nodes, elements, sections, and - forces for the beam are defined in the analysis model - - Note, nodes and elements will both start at 0 instead of 1. - - Parameters - ---------- - beam : planesections Beam2D - The beam whose data is being recorded. - recorder : planesections Recorder - The recorder to use for the output beam. - geomTransform: str, optional - The OpenSees Geometry transform to use. Can be "Linear" or "PDelta" - clearOld : bool, optional - A flag that can be used to turn on or off clearing the old analysis - when the beam is created. - There are some very niche cases where users may want to have mutiple - beams at once in the OpenSees model. - However, this should remain true for nearly all analyses. - Do not turn on unless you know what you're doing. - - """ - - def __init__(self, beam3D:bb.Beam, recorder = OutputRecorderOpenSees, - geomTransform = 'Linear', clearOld = True): - - self.beam = beam3D - self._checkBeam(beam3D) - - self.Nnode = beam3D.Nnodes - self.Nele = self.Nnode - 1 - self.recorder = OutputRecorderOpenSees - self.geomTransform = geomTransform - self.clearOld = clearOld - - def _checkBeam(self, beam3D): - if not beam3D._dimension: - raise Exception("The beam has no dimension, something terible has happened.") - if beam3D._dimension != '3D': - raise Exception(f"The beam has dimension of type {beam3D._dimension}, it should have type '3D'") - - def initModel(self, clearOld = True): - """ - Initializes the model. - - Parameters - ---------- - clearOld : bool, optional - A flag that can be used to turn on or off clearing the old analysis - when the beam is created. - There are some very niche cases where users may want to have mutiple - beams at once in the OpenSees model. - However, this should remain true for nearly all analyses. - Do not turn on unless you know what you're doing. - """ - - if clearOld: - op.wipe() - op.model('Basic' , '-ndm', 3) - # see https://opensees.berkeley.edu/wiki/index.php/PDelta_Transformation - op.geomTransf(self.geomTransform, 1, *[0, 0, 1]) - - - def buildNodes(self): - """ - Adds each node in the beam to the OpenSeesPy model, and assigns - that node a fixity. - """ - - for node in self.beam.nodes: - op.node(int(node.ID), float(node.x), 0., 0.) - - # OpenSees is very finicky with these inputs, int them for saftey. - f1, f2, f3, f4, f5, f6 = node.fixity.fixityValues - op.fix(node.ID, int(f1), int(f2), int(f3), int(f4), int(f5), int(f6)) - - def buildEulerBeams(self): - """ - Creates an elastic Euler beam between each node in the model. - """ - beam = self.beam - matPropreties = beam.getMaterialPropreties() - for ii in range(self.Nele): - ID = ii + 1 - eleID = int(ID) - Ni = int(ID) - Nj = int(ID + 1) - # element('elasticBeamColumn', eleTag, *eleNodes, Area, E_mod, G_mod, Jxx, Iy, Iz, transfTag, <'-mass', mass>, <'-cMass'>) - op.element(beam.EleType, eleID , Ni, Nj , *matPropreties, 1) - - def buildPointLoads(self): - """ - Applies point loads to the appropriate nodes in the model. - """ - op.timeSeries('Linear',1) - op.pattern('Plain', 1, 1) - for load in self.beam.pointLoads: - op.load(int(load.nodeID), *load.P) - - - def buildAnalysisPropreties(self): - """ - Typical openSeesPy propreties that should work for any linear beam. - A linear algorithm is used because there is no nonlienarity in the beam. - """ - # op.constraints("Transformation") - op.constraints("Lagrange") - op.numberer("Plain") - op.system('BandGeneral') - op.test('NormDispIncr', 1.*10**-8, 40, 0 , 2) - # op.algorithm('Newton') - op.algorithm('Linear') - op.integrator('LoadControl',1.) - op.analysis('Static') - - def analyze(self): - """ - Analyzes the model once and records outputs. - """ - ok = op.analyze(1) - op.reactions() - - - def buildEleLoads(self): - """ - Applies element loads to the appropriate elements in the model. - """ - op.timeSeries('Linear',2) - op.pattern('Plain', 2, 2) - - for eleload in self.beam.eleLoads: - N1 = self.beam._findNode(eleload.x1) + 1 - N2 = self.beam._findNode(eleload.x2) + 1 - load = eleload.P - - for ii in range(N1, N2): - # eleLoad('-ele', *eleTags, '-range', eleTag1, eleTag2, '-type', '-beamUniform', Wy, , Wx=0.0, '-beamPoint',Py,,xL,Px=0.0,'-beamThermal',*tempPts) - op.eleLoad('-ele', int(ii), '-type', '-beamUniform', load[1], load[2], load[0]) - - def runAnalysis(self, recordOutput = True): - """ - Makes and analyzes the beam in OpenSees. - - Returns - ------- - None. - - """ - - self.initModel(self.clearOld) - self.buildNodes() - self.buildEulerBeams() - self.buildPointLoads() - self.buildEleLoads() - self.buildAnalysisPropreties() - self.analyze() - - if recordOutput == True: - self.recorder(self.beam) \ No newline at end of file diff --git a/src/builder.py b/src/builder.py index 71abe94..1a5eba8 100644 --- a/src/builder.py +++ b/src/builder.py @@ -679,15 +679,20 @@ def addPointLoad(self, x:float, pointLoad:list, label:str = '', A vertical load of 10 is applied in beam units. New Load 2: [0., 0., 13] - New node 2: A moment of 13 is applied in beam units. + A moment of 13 is applied in beam units. label : str, optional The label of the input node. labels are displayed in the plots. The default is ''. """ + # Catch incorrectly given types. if hasattr(pointLoad, '__iter__') == False: raise Exception('Point load vector must be a list or Numpy array.') + + # Converty to np array. We do vector operations on data downstream + if isinstance(pointLoad, list): + pointLoad = np.array(pointLoad) loadID = len(self.pointLoads) + 1 diff --git a/src/postprocess/__init__.py b/src/postprocess/__init__.py index aa38884..d34bea5 100644 --- a/src/postprocess/__init__.py +++ b/src/postprocess/__init__.py @@ -1,8 +1,4 @@ -from .plot import (getDisp, getVertDisp, getMaxVertDisp, - plotInternalForce, plotShear, plotMoment, +from .parse import getDisp, getVertDisp, getMaxVertDisp +from .plot import (plotInternalForce, plotShear, plotMoment, plotDisp, plotVertDisp, plotRotation, getInternalForces2D, plotMoment2D, plotShear2D) - - - - diff --git a/src/postprocess/parse.py b/src/postprocess/parse.py new file mode 100644 index 0000000..a03dea4 --- /dev/null +++ b/src/postprocess/parse.py @@ -0,0 +1,98 @@ +""" +Common post-processing functions for plotting. +""" + +import numpy as np +from planesections import Beam + + + + +def getDisp(beam: Beam, ind: int): + """ + Gets the beam displacement along the axis specified for the index. + + Parameters + ---------- + beam : Beam + The beam to read displacement from. The beam must be analyzed to get + data. + ind : int + The index of the axis to read from. Can have value 0: horizontal + displacement + 1: vertical displacement + 2: rotation. + + Returns + ------- + disp : numpy array + The displacement at each x coordinant. + xcoords : numpy array + The x coordinants. + """ + + xcoords = np.zeros(beam.Nnodes) + disp = np.zeros(beam.Nnodes) + for ii, node in enumerate(beam.nodes): + xcoords[ii] = node.x + disp[ii] = node.disps[ind] + return disp, xcoords + + + +def getVertDisp(beam: Beam): + """ + Gets the beam vertical displacement for the beam + + Parameters + ---------- + beam : Beam + The beam to read displacement from. The beam must be analyzed to get + data. + + Returns + ------- + disp : numpy array + The displacement at each x coordinant. + xcoords : numpy array + The x coordinants. + """ + return getDisp(beam, 1) + + +def getMaxVertDisp(beam: Beam): + """ + Gets the absolute value of beam vertical displacement and it's location. + + Parameters + ---------- + beam : Beam + The beam to read displacement from. The beam must be analyzed to get + data. + + Returns + ------- + dispMax : float + The displacement at each x coordinant. + xcoords : numpy array + The x coordinants. + """ + disp, x = getVertDisp(beam) + dispAbs = np.abs(disp) + ind = np.argmax(dispAbs) + return disp[ind], x[ind] + + +def _getForceValues(beam, index): + Nnodes = len(beam.nodes) + xcoords = np.zeros(Nnodes*2) + force = np.zeros(Nnodes*2) + labels = [None]*Nnodes + for ii, node in enumerate(beam.nodes): + ind1 = 2*ii + ind2 = ind1 + 2 + xcoords[ind1:ind2] = node.x + force[ind1:ind2] = node.getInternalForces(index) + labels[ii] = node.label + + return xcoords, force, labels \ No newline at end of file diff --git a/src/postprocess/plot.py b/src/postprocess/plot.py index c6bf109..12c90a8 100644 --- a/src/postprocess/plot.py +++ b/src/postprocess/plot.py @@ -1,108 +1,14 @@ import numpy as np import matplotlib.pyplot as plt from planesections import Node2D, Beam +import textalloc as ta +from dataclasses import dataclass -# ============================================================================= -# -# ============================================================================= - - -def getDisp(beam: Beam, ind: int): - """ - Gets the beam displacement along the axis specified for the index. - - Parameters - ---------- - beam : Beam - The beam to read displacement from. The beam must be analyzed to get data. - ind : int - The index of the axis to read from. Can have value 0: horizontal displacement - 1: vertical displacement - 2: rotation. - - Returns - ------- - disp : numpy array - The displacement at each x coordinant. - xcoords : numpy array - The x coordinants. - """ - - xcoords = np.zeros(beam.Nnodes) - disp = np.zeros(beam.Nnodes) - for ii, node in enumerate(beam.nodes): - xcoords[ii] = node.x - disp[ii] = node.disps[ind] - return disp, xcoords - - - -def getVertDisp(beam: Beam): - """ - Gets the beam vertical displacement for the beam - - Parameters - ---------- - beam : Beam - The beam to read displacement from. The beam must be analyzed to get data. - - Returns - ------- - disp : numpy array - The displacement at each x coordinant. - xcoords : numpy array - The x coordinants. - """ - return getDisp(beam, 1) - - -def getMaxVertDisp(beam: Beam): - """ - Gets the absolute value of beam vertical displacement and it's location. - - Parameters - ---------- - beam : Beam - The beam to read displacement from. The beam must be analyzed to get data. - - Returns - ------- - dispMax : float - The displacement at each x coordinant. - xcoords : numpy array - The x coordinants. - """ - disp, x = getVertDisp(beam) - dispAbs = np.abs(disp) - ind = np.argmax(dispAbs) - return disp[ind], x[ind] +from .parse import _getForceValues +from .poi import findAllPOI, removeFalsePOI - - - - - - -# def getMaxBeamDisp(beam: Beam, ind: int): -# """ -# Gets the beam displacement along the axis specified for the index. - -# Parameters -# ---------- -# beam : Beam -# DESCRIPTION. -# ind : TYPE -# DESCRIPTION. - -# Returns -# ------- -# None. - -# """ - - # ============================================================================= # Plotting # ============================================================================= @@ -191,7 +97,22 @@ def plotMoment(beam:Beam, scale:float=-1, yunit = 'Nm', **kwargs): Turns on or off the grid. labelPlot : bool, optional Turns on or off the plot labels. + labelPOI : bool, optional + Turns on or off the plot point of interst labels. + POIOptions : dict, optional + The options to use for the POI labels. There are several flags that + have values true or false which can be used to turn on or off certain + points of interst. Notably: + + showLabels: this flag turns on or off + + showDiscontinutiy: this flag turns on labeling points of discontinuity + + showMax: this turns on or off + The default is None, which sets all flags to true + + Returns ------- fig : matplotlib fig @@ -245,7 +166,7 @@ def plotShear(beam:Beam, scale:float=1, **kwargs): Kwarg Parameters ---------- - These are possible Kwarg parameters that will be passed to plotInternalForce2D. + These are possible Kwarg parameters that will be passed to plotInternalForce. xunit : str, optional The xunit for the plot labels. The default is m. @@ -257,6 +178,21 @@ def plotShear(beam:Beam, scale:float=1, **kwargs): Turns on or off the grid. labelPlot : bool, optional Turns on or off the plot labels. + labelPOI : bool, optional + Turns on or off the plot point of interst labels. + POIOptions : dict, optional + The options to use for the POI labels. There are several flags that + have values true or false which can be used to turn on or off certain + points of interst. Notably: + + showLabels: this flag turns on or off + + showDiscontinutiy: this flag turns on labeling points of discontinuity + + showMax: this turns on or off + + The default is None, which sets all flags to true + Returns ------- @@ -265,30 +201,22 @@ def plotShear(beam:Beam, scale:float=1, **kwargs): ax : matplotlib ax line : matplotlib line - the plotted line. + The plotted line. """ return plotInternalForce(beam, 1, scale, **kwargs) - -def _getForceValues(beam, index): - Nnodes = len(beam.nodes) - xcoords = np.zeros(Nnodes*2) - force = np.zeros(Nnodes*2) - labels = [None]*Nnodes - for ii, node in enumerate(beam.nodes): - ind1 = 2*ii - ind2 = ind1 + 2 - xcoords[ind1:ind2] = node.x - force[ind1:ind2] = node.getInternalForces(index) - labels[ii] = node.label - - return xcoords, force, labels +@dataclass +class _NodeOutputs: + xcoords:list[float] + force:list[float] + labels:list[str] def plotInternalForce(beam:Beam, index:int, scale:float, xunit= 'm', yunit = 'N', - showAxis = True, showGrid = False, labelPlot = True): + showAxis = True, showGrid = False, labelPlot = True, + labelPOI = False, POIOptions = None): """ Plots the internal forces within a beam. Two values are given for each point, one at the left side, one at the right @@ -322,7 +250,20 @@ def plotInternalForce(beam:Beam, index:int, scale:float, xunit= 'm', yunit = 'N' Turns on or off the grid. labelPlot : bool, optional Turns on or off the plot labels. - + labelPOI : bool, optional + Turns on or off the plot point of interst labels. + POIOptions : dict, optional + The options to use for the POI labels. There are several flags that + have values true or false which can be used to turn on or off certain + points of interst. Notably: + showLabels: this flag turns on or off POI label on points with user labels. + + showDiscontinutiy: this flag turns on POI labels on points of discontinuity. + + showMax: this turns on or off POI at maximum points. + + The default is None, which sets all flags to true + Returns ------- fig : matplotlib fig @@ -334,24 +275,30 @@ def plotInternalForce(beam:Beam, index:int, scale:float, xunit= 'm', yunit = 'N' """ xcoords, force, labels = _getForceValues(beam, index) - + forceScaled = force*scale + fig, ax = _initOutputFig(showAxis, showGrid) _plotAxis(ax, xcoords, xunit, yunit) if labelPlot: - _plotLabels(ax, xcoords[::2], force*scale, labels) + _plotLabels(ax, xcoords[::2], forceScaled, labels) + line = plt.plot(xcoords, forceScaled) - line = plt.plot(xcoords, force*scale) + if labelPOI: + shear = None + if index ==2: + _, shear, _ = _getForceValues(beam, index) + candidatePOI = findAllPOI(xcoords, forceScaled, labels, shear, POIOptions) + filteredPoiInd = removeFalsePOI(candidatePOI, force) + plotPOI(fig, ax, xcoords, forceScaled, labels, filteredPoiInd) return fig, ax, line def plotInternalForce2D(beam:Beam, index:int, scale:float, xunit= 'm', yunit = 'N', showAxis = True, showGrid = False, labelPlot = True): - print('All 2D plot functions are depricated and will return an error in future releases. \ + raise Exception('All 2D plot functions are depricated and will be removed in future releases. \ Use the generic plot functions instead, i.e. plotInternalForce instead of plotInternalForce2D' ) - return plotInternalForce(beam, index, scale, 'm', 'N', - True, False, True) def plotVertDisp(beam:Beam, scale=1000, yunit = 'mm', **kwargs): @@ -432,12 +379,9 @@ def plotRotation(beam:Beam, scale=1000, yunit = 'mRad', **kwargs): return plotDisp(beam, ind, scale, yunit= yunit, **kwargs) - - - - def plotDisp(beam:Beam, index=1, scale=1000, xunit= 'm', yunit = 'mm', - showAxis = True, showGrid = False, labelPlot = True): + showAxis = True, showGrid = False, labelPlot = True, + labelPOI = False, POIOptions = None): """ Plots the displacement of one dimension of the beam. Each node will contain the relevant dispancement information. Analysis must be run on the beam prior to @@ -485,15 +429,19 @@ def plotDisp(beam:Beam, index=1, scale=1000, xunit= 'm', yunit = 'mm', fig, ax = _initOutputFig(showAxis, showGrid) _plotAxis(ax, xcoords, xunit, yunit, 'Displacement') - + disp = disp*scale + if labelPlot: - _plotLabels(ax, xcoords, disp*scale, labels) + _plotLabels(ax, xcoords, disp, labels) - line = plt.plot(xcoords, disp*scale) - - return fig, ax, line + line = plt.plot(xcoords, disp) + if labelPOI: + candidatePOI = findAllPOI(xcoords, disp, labels, POIOptions = POIOptions) + filteredPoiInd = removeFalsePOI(candidatePOI, disp) + plotPOI(fig, ax, xcoords, disp, labels, filteredPoiInd) + return fig, ax, line def plotDisp2D(beam:Beam, index=1, scale=1000, xunit= 'm', yunit = 'mm', showAxis = True, showGrid = False, labelPlot = True): @@ -509,7 +457,6 @@ def plotVertDisp2D(beam:Beam, scale=1000, yunit = 'mm', **kwargs): return plotDisp2D(beam, 1, scale, yunit= yunit, **kwargs) - def plotRotation2D(beam:Beam, scale=1000, yunit = 'mRad', **kwargs): """ Depricated, instead use plotRotation @@ -520,3 +467,31 @@ def plotRotation2D(beam:Beam, scale=1000, yunit = 'mRad', **kwargs): return plotDisp2D(beam, ind, scale, yunit= yunit, **kwargs) + +def plotPOI(fig, ax, xcoords, force, labels, filteredPoiInd): + + labelX = [] + labelY = [] + labelText = [] + for ind in filteredPoiInd: + labelName = labels[int((ind+1)/2)] + x0 = xcoords[ind] + y0 = force[ind] + xOut = round(x0,2) + yOut = round(y0,1) + base = '' + if labelName: + base = ' Point ' + labelName + ' \n' + textX = f' x = {xOut} \n' + textY = f' y = {yOut} ' + text = base + textX + textY + labelX.append(x0) + labelY.append(y0) + labelText.append(text) + + ta.allocate_text(fig, ax, labelX, labelY, labelText, textsize=8, + x_scatter=labelX, y_scatter=labelY, + x_lines=[xcoords], y_lines=[force], + linecolor='grey', linewidth=0.5, + max_distance = 0.5, min_distance = 0.02, + seed =1,) \ No newline at end of file diff --git a/src/postprocess/poi.py b/src/postprocess/poi.py index 2fb07ad..eff50b4 100644 --- a/src/postprocess/poi.py +++ b/src/postprocess/poi.py @@ -1,4 +1,297 @@ import numpy as np -import matplotlib.pyplot as plt -from planesections import Node2D, Beam +from .parse import _getForceValues +def _get_discontinuity_inds(force, tol = 10e-6): + + force2D = force.reshape(-1,2) + indDisTmp = np.where(tol < np.abs(np.diff(force2D)))[0] + indDisTmp = np.concatenate((indDisTmp*2, indDisTmp*2 + 1)) + # indDis = indDisTmp[1:-1] + return indDisTmp + +def _valsInPercentTol(y1, y2, tolPercent = 10e-4): + if abs((y1/y2) - 1) < tolPercent: + return True + return False + +def valRelativelyNearZero(val, referenceMax, tolPercent = 10e-9): + """ + Remove The final point if it is close to zero. + We do some gymnastics here because we don't know the scale of the plot. + """ + return (abs(val / referenceMax) < tolPercent) + + + +def _getLabelInds(xcoords:list, force:list, labels:list): + nonNoneLabels = np.array(labels) != None + nonEmptyLabels = np.array(labels) != '' + + indLabelTmp = np.where(nonNoneLabels*nonEmptyLabels)[0] + + indLabelTmp = indLabelTmp*2 + indLabelTmp = np.array([indLabelTmp, indLabelTmp + 1]).T + + # If the label index is at a discontinous point, get both sides of the point. + indLabel = [] + tol = 10e-6 + for ind1, ind2 in indLabelTmp: + xcoordEqual = (xcoords[ind1] == xcoords[ind2]) + ycoordUnequal = tol < abs(force[ind1]/force[ind2] - 1) + if xcoordEqual and ycoordUnequal: + indLabel += [ind1, ind2] + else: + indLabel.append(ind1) + return indLabel + + +def _getLabelIndsDisp(xcoords:list, force:list, labels:list): + nonNoneLabels = np.array(labels) != None + nonEmptyLabels = np.array(labels) != '' + + indLabelTmp = np.where(nonNoneLabels*nonEmptyLabels)[0] + indLabel = np.array(indLabelTmp,dtype =int) + + return indLabel + +def _processPOIDictInputs(poiOptions): + + if 'showLabels' not in poiOptions: + poiOptions['showLabels'] = True + + if 'showDiscontinutiy' not in poiOptions: + poiOptions['showDiscontinutiy'] = True + + if 'showMax' not in poiOptions: + poiOptions['showMax'] = True + + + return poiOptions + + + +def _getPOIIndsForce(xcoords, ycoords, labels, ycordsSecondary, poiOptions): + + # Get Label Inds + indLabel = [] + if poiOptions['showLabels'] == True: + indLabel = _getLabelInds(xcoords, ycoords, labels) + + # Get Discontinuity Indicies + indDis = [] + if poiOptions['showDiscontinutiy'] == True: + indDis = _get_discontinuity_inds(ycoords) + + # If it's a moment plot, catch points of shear discontinuty. + indDisShear = [] + if ycordsSecondary is not None: + indDisShear = _get_discontinuity_inds(ycordsSecondary) + + # Get Max/min Indicies + maxMinInds = [] + if poiOptions['showMax'] == True: + indMax = np.argmax(ycoords) + indMin = np.argmin(ycoords) + maxMinInds = [indMax, indMin] + + return np.concatenate((indDis, indDisShear, np.array(indLabel, dtype=int), maxMinInds)) + + + +def _getPOIIndsDisp(xcoords, ycoords, labels, ycordsSecondary, poiOptions): + + # Get Label Inds + indLabel = [] + if poiOptions['showLabels'] == True: + indLabel = _getLabelIndsDisp(xcoords, ycoords, labels) + + # Get Max/min Indicies + maxMinInds = [] + if poiOptions['showMax'] == True: + indMax = np.argmax(ycoords) + indMin = np.argmin(ycoords) + maxMinInds = [indMax, indMin] + + return np.concatenate((np.array(indLabel, dtype=int), maxMinInds)) + + + +def findBeamForcePOI(beam, index, POIOptions): + """ + Gets the Force indicies of all the points of interest for the input beam. + + Parameters + ---------- + beam : ndarray + The input beam. + index : int + The type of response to plot, can have value + + In 2D has values: + [0:Fx, 1: Fy, 2: M] + + In 3D has values: + [0:Fx, 1: Fy, 2: Fz, 3: Mx, 4: Mx, 5: Mz] + + + POIOptions : dict, optional + The options to use for the POI labels. There are several flags that + have values true or false which can be used to turn on or off certain + points of interst. Notably: + + showLabels: this flag turns on or off + + showDiscontinutiy: this flag turns on labeling points of discontinuity + + showMax: this turns on or off + + The default is None, which sets all flags to true + + Returns + ------- + poiInd : list[int] + The list of indexes for each point of interest. + + """ + + xcoords, force, labels = _getForceValues(beam, index) + shear = None + if index ==2: + _, shear, _ = _getForceValues(beam, index) + + candidatePOI = findAllPOI(xcoords, force, labels, shear, POIOptions) + return removeFalsePOI(candidatePOI, force) + + + +def findAllPOI(xcoords, ycoords, labels, ysecondary=None, POIOptions:dict = None): + """ + Gets the indicies of all the points of interest for the input beam. + + Parameters + ---------- + xcoords : ndarray + The input x coordinant array. + ycoords : ndarray + The input y coordinant array. + labels : list + The input list of labels for each node. + ysecondary : TYPE, optional + A secondary array that can be used to find discontinous points. + The default is None. + POIOptions : dict, optional + The options to use for the POI labels. There are several flags that + have values true or false which can be used to turn on or off certain + points of interst. Notably: + showLabels: this flag turns on or off + showDiscontinutiy: this flag turns on labeling points of discontinuity + showMax: this turns on or off + + The default is None, which sets all flags to true + + Returns + ------- + poiInd : list[int] + The list of indexes for each point of interest. + + """ + + # Prevent the options from having no value. + if POIOptions == None: + POIOptions = {} + + # add default values if they aren't included + POIOptions = _processPOIDictInputs(POIOptions) + + # Check if the + if len(ycoords) == len(labels): + union = _getPOIIndsDisp(xcoords, ycoords, labels, ysecondary, POIOptions) + else: + union = _getPOIIndsForce(xcoords, ycoords, labels, ysecondary, POIOptions) + + poiInd = list(union) + poiInd = [int(ind) for ind in poiInd] + return poiInd + +# ============================================================================= +# +# ============================================================================= + + +def _findCloseDiscontinousPoints(force, ind, poiInd): + """ + This will remove discontinous points picked up in the shear force diagram. + + We only filter out the left most point, otherwise we would end up filtering + both points out. + + """ + # Combine discontinuous points that are within a tolerance of eachother + leftAdjecentPointExists = False + for indOther in poiInd: + if (indOther - 1 == ind ): + leftAdjecentPointExists = True + break + + # don't add the point if it is next to another point and close in value + if leftAdjecentPointExists and _valsInPercentTol(force[ind], force[indOther]): + return True + else: + return False + +def removeFalsePOI(candidatePOI, force): + """ + Do a final pruning of points + + + Parameters + ---------- + beam : TYPE + DESCRIPTION. + forceInd : TYPE + DESCRIPTION. + candidatePOI : TYPE + DESCRIPTION. + + Returns + ------- + None. + + """ + + # get the forces again + absMax = np.max(np.abs(force)) + + # removeThe end points + end = len(force) - 1 + start = 0 + if end in candidatePOI: + candidatePOI.remove(end) + if start in candidatePOI: + candidatePOI.remove(start) + + # make a list of the second from last points. + # If these are close to zero we will remove them later. + candidateEndPoints = [1, len(force) - 2] + + filteredPoiInd = [] + poiInd = [int(ind) for ind in candidatePOI] + for ind in poiInd: + ind = int(ind) + + if _findCloseDiscontinousPoints(force, ind, poiInd): + continue + + # Values + if ind in candidateEndPoints and valRelativelyNearZero(force[ind], absMax): + continue + + filteredPoiInd.append(ind) + filteredPoiInd = set(filteredPoiInd) + return filteredPoiInd + + + +# ============================================================================= +# +# ============================================================================= diff --git a/src/requirements.txt b/src/requirements.txt index bc6cec5..cc614ee 100644 --- a/src/requirements.txt +++ b/src/requirements.txt @@ -1,3 +1,5 @@ matplotlib numpy openseespy +pynite +textalloc>=0.0.6