From 0bf0543a66de2252b8bf634fc5ffb66d64af5593 Mon Sep 17 00:00:00 2001 From: claudio perez <50180406+claudioperez@users.noreply.github.com> Date: Sat, 31 Aug 2024 21:03:43 -0700 Subject: [PATCH] cmp - merging repos --- {ops => SRC}/opensees/README.md | 0 {ops => SRC}/opensees/__init__.py | 0 {ops => SRC}/opensees/__main__.py | 0 SRC/opensees/elements.py | 953 ------- {ops => SRC}/opensees/emit/__init__.py | 0 {ops => SRC}/opensees/emit/__main__.py | 0 {ops => SRC}/opensees/emit/apidoc.py | 0 {ops => SRC}/opensees/emit/emit_json.py | 0 {ops => SRC}/opensees/emit/emitter.py | 0 {ops => SRC}/opensees/emit/fedeas.py | 0 {ops => SRC}/opensees/emit/mesh.py | 0 {ops => SRC}/opensees/emit/opensees.py | 0 {ops => SRC}/opensees/emit/registry.py | 0 {ops => SRC}/opensees/emit/writer.py | 0 {ops => SRC}/opensees/library/__init__.py | 0 {ops => SRC}/opensees/library/ast.py | 0 {ops => SRC}/opensees/library/material.py | 0 {ops => SRC}/opensees/library/model.py | 0 {ops => SRC}/opensees/library/obj.py | 0 {ops => SRC}/opensees/library/pattern.py | 0 SRC/opensees/material.py | 1 - SRC/opensees/matrices.py | 2333 ----------------- {ops => SRC}/opensees/repl/__init__.py | 0 {ops => SRC}/opensees/repl/cmdshell.py | 0 {ops => SRC}/opensees/repl/ptkshell.py | 0 {ops => SRC}/opensees/repl/unix.py | 0 SRC/opensees/scripts.py | 53 - {ops => SRC}/opensees/section/__init__.py | 0 .../opensees/section/aisc_imperial.py | 0 {ops => SRC}/opensees/section/analysis.py | 0 {ops => SRC}/opensees/section/patch.py | 0 {ops => SRC}/opensees/section/section.py | 0 {ops => SRC}/opensees/section/torsion.py | 0 SRC/opensees/state.py | 4 - {ops => SRC}/opensees/units/__init__.py | 0 {ops => SRC}/opensees/units/__main__.py | 0 {ops => SRC}/opensees/units/cgs.py | 0 {ops => SRC}/opensees/units/english.py | 0 {ops => SRC}/opensees/units/fas.py | 0 {ops => SRC}/opensees/units/fks.py | 0 {ops => SRC}/opensees/units/fps.py | 0 {ops => SRC}/opensees/units/ias.py | 0 {ops => SRC}/opensees/units/iks.py | 0 {ops => SRC}/opensees/units/ips.py | 0 {ops => SRC}/opensees/units/mks.py | 0 {ops => SRC}/opensees/units/si.py | 0 {ops => SRC}/opensees/units/units.py | 0 SRC/opensees/validation.py | 11 - COPYRIGHT => about/COPYRIGHT | 0 ops/opensees/analysis.py | 134 - ops/opensees/openseespy.py | 790 ------ ops/opensees/recorder.py | 177 -- ops/opensees/series.py | 34 - ops/opensees/tcl.py | 390 --- .clang-format => tools/.clang-format | 0 tools/SCRIPTS/strain.tcl | 389 +++ tools/SCRIPTS/toOpenSeesPy.py | 116 + tools/SCRIPTS/uniaxialMaterials.tcl | 27 + tools/apidoc.lua | 30 + tools/build-conda-package.sh | 52 + tools/doc.py | 173 ++ tools/docopt_c.py | 701 +++++ tools/get_wiki_img.sh | 7 + tools/install.sh | 10 + tools/pandoc_plantuml_filter.py | 51 + tools/tcl2py.sh | 159 ++ 66 files changed, 1715 insertions(+), 4880 deletions(-) rename {ops => SRC}/opensees/README.md (100%) rename {ops => SRC}/opensees/__init__.py (100%) rename {ops => SRC}/opensees/__main__.py (100%) delete mode 100644 SRC/opensees/elements.py rename {ops => SRC}/opensees/emit/__init__.py (100%) rename {ops => SRC}/opensees/emit/__main__.py (100%) rename {ops => SRC}/opensees/emit/apidoc.py (100%) rename {ops => SRC}/opensees/emit/emit_json.py (100%) rename {ops => SRC}/opensees/emit/emitter.py (100%) rename {ops => SRC}/opensees/emit/fedeas.py (100%) rename {ops => SRC}/opensees/emit/mesh.py (100%) rename {ops => SRC}/opensees/emit/opensees.py (100%) rename {ops => SRC}/opensees/emit/registry.py (100%) rename {ops => SRC}/opensees/emit/writer.py (100%) rename {ops => SRC}/opensees/library/__init__.py (100%) rename {ops => SRC}/opensees/library/ast.py (100%) rename {ops => SRC}/opensees/library/material.py (100%) rename {ops => SRC}/opensees/library/model.py (100%) rename {ops => SRC}/opensees/library/obj.py (100%) rename {ops => SRC}/opensees/library/pattern.py (100%) delete mode 100644 SRC/opensees/material.py delete mode 100644 SRC/opensees/matrices.py rename {ops => SRC}/opensees/repl/__init__.py (100%) rename {ops => SRC}/opensees/repl/cmdshell.py (100%) rename {ops => SRC}/opensees/repl/ptkshell.py (100%) rename {ops => SRC}/opensees/repl/unix.py (100%) delete mode 100644 SRC/opensees/scripts.py rename {ops => SRC}/opensees/section/__init__.py (100%) rename {ops => SRC}/opensees/section/aisc_imperial.py (100%) rename {ops => SRC}/opensees/section/analysis.py (100%) rename {ops => SRC}/opensees/section/patch.py (100%) rename {ops => SRC}/opensees/section/section.py (100%) rename {ops => SRC}/opensees/section/torsion.py (100%) delete mode 100644 SRC/opensees/state.py rename {ops => SRC}/opensees/units/__init__.py (100%) rename {ops => SRC}/opensees/units/__main__.py (100%) rename {ops => SRC}/opensees/units/cgs.py (100%) rename {ops => SRC}/opensees/units/english.py (100%) rename {ops => SRC}/opensees/units/fas.py (100%) rename {ops => SRC}/opensees/units/fks.py (100%) rename {ops => SRC}/opensees/units/fps.py (100%) rename {ops => SRC}/opensees/units/ias.py (100%) rename {ops => SRC}/opensees/units/iks.py (100%) rename {ops => SRC}/opensees/units/ips.py (100%) rename {ops => SRC}/opensees/units/mks.py (100%) rename {ops => SRC}/opensees/units/si.py (100%) rename {ops => SRC}/opensees/units/units.py (100%) delete mode 100644 SRC/opensees/validation.py rename COPYRIGHT => about/COPYRIGHT (100%) delete mode 100644 ops/opensees/analysis.py delete mode 100644 ops/opensees/openseespy.py delete mode 100644 ops/opensees/recorder.py delete mode 100644 ops/opensees/series.py delete mode 100644 ops/opensees/tcl.py rename .clang-format => tools/.clang-format (100%) create mode 100755 tools/SCRIPTS/strain.tcl create mode 100644 tools/SCRIPTS/toOpenSeesPy.py create mode 100644 tools/SCRIPTS/uniaxialMaterials.tcl create mode 100644 tools/apidoc.lua create mode 100644 tools/build-conda-package.sh create mode 100644 tools/doc.py create mode 100644 tools/docopt_c.py create mode 100755 tools/get_wiki_img.sh create mode 100644 tools/install.sh create mode 100755 tools/pandoc_plantuml_filter.py create mode 100755 tools/tcl2py.sh diff --git a/ops/opensees/README.md b/SRC/opensees/README.md similarity index 100% rename from ops/opensees/README.md rename to SRC/opensees/README.md diff --git a/ops/opensees/__init__.py b/SRC/opensees/__init__.py similarity index 100% rename from ops/opensees/__init__.py rename to SRC/opensees/__init__.py diff --git a/ops/opensees/__main__.py b/SRC/opensees/__main__.py similarity index 100% rename from ops/opensees/__main__.py rename to SRC/opensees/__main__.py diff --git a/SRC/opensees/elements.py b/SRC/opensees/elements.py deleted file mode 100644 index 608e80a1ec..0000000000 --- a/SRC/opensees/elements.py +++ /dev/null @@ -1,953 +0,0 @@ -import numpy as np -from abc import abstractmethod - -from numpy.polynomial import Polynomial -from scipy.integrate import quad -import scipy.integrate - -from .matrices import Structural_Matrix, Structural_Vector - - -class IntForce: - def __init__(self, elem, number): - self.number = number - self.elem = elem - self.rel = False - self.redundant = False - self.plastic_event = None - - @property - def tag(self): - return self.elem.tag + '_' + str(self.number) - - def __str__(self): - return self.tag - - def __repr__(self): - return self.tag - - -class BasicLink(): - """Class implementing general geometric element methods""" - - def __init__(self, ndf, ndm, nodes): - self.nodes = nodes - self.ndf: int = ndf - self.ndm: int = ndm - self.nen = len(nodes) - - @property - def dofs(self): - eid = np.array([]) - for node in self.nodes: - eid = np.append(eid, node.dofs[0:self.ndf]) - return eid - - @property - def L(self): - xyzi = self.nodes[0].xyz - xyzj = self.nodes[1].xyz - L = np.linalg.norm(xyzi-xyzj) - return L - - @property - def L0(self): - xyzi = self.nodes[0].xyz0 - xyzj = self.nodes[1].xyz0 - L = np.linalg.norm(xyzi-xyzj) - return L - - @property - def Li(self): - n1 = self.nodes[0] - n2 = self.nodes[1] - xyzi = np.array([n1.xi, n1.yi, n1.zi]) - xyzj = np.array([n2.xi, n2.yi, n2.zi]) - L = np.linalg.norm(xyzi-xyzj) - return L - - - @property - def Dx(self): - return self.nodes[1].x - self.nodes[0].x - - @property - def Dy(self): - return self.nodes[1].y - self.nodes[0].y - - @property - def Dz(self): - return self.nodes[1].z - self.nodes[0].z - - @property - def sn(self): - L = self.L - sn = (self.nodes[1].y - self.nodes[0].y)/L - return sn - - @property - def cs(self): - L = self.L - cs = (self.nodes[1].x - self.nodes[0].x)/L - return cs - - @property - def cz(self): - L = self.L - cz = (self.nodes[1].z - self.nodes[0].z)/L - return cz - - def Rx_matrix(self): - """Rotation about x - - https://en.wikipedia.org/wiki/Rotation_matrix - """ - cs = self.cs - sn = self.sn - Rx = np.array([ - [1.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [0.0, cs, -sn, 0.0, 0.0, 0.0], - [0.0, sn, cs, 0.0, 0.0, 0.0], - [0.0, 0.0, 0.0, 1.0, 0.0, 0.0], - [0.0, 0.0, 0.0, 0.0, cs, -sn], - [0.0, 0.0, 0.0, 0.0, sn, cs], - ]) - - return 0 - - def Ry_matrix(self): - "Rotation about z" - cs = self.cs - sn = self.sn - Ry = np.array([ - [ cs, 0.0, sn], - [0.0, 1.0, 0.0], - [-sn, 0.0, cs], - ]) - - return 0 - - def Rz_matrix(self): - "Rotation about z" - cs = self.cs - sn = self.sn - Rz = np.array([ - [ cs, -sn, 0.0, 0.0, 0.0, 0.0], - [ sn, cs, 0.0, 0.0, 0.0, 0.0], - [0.0, 0.0, 1.0, 0.0, 0.0, 0.0], - [0.0, 0.0, 0.0, cs, -sn, 0.0], - [0.0, 0.0, 0.0, sn, cs, 0.0], - [0.0, 0.0, 0.0, 0.0, 0.0, 1.0]]) - - return Rz - - -class Element(BasicLink): - """Element parent class""" - - def __init__(self, ndf, ndm, force_dict, nodes): - super().__init__( ndf, ndm, nodes) - - self.history = {} - self.current = {} - - nq = len(force_dict) - self.rel = {str(i): False for i in range(1, nq+1)} - self.red = {str(i): False for i in range(1, nq+1)} - self.q0 = {str(i): 0. for i in range(1, nq+1)} - self.v0 = {str(i): 0. for i in range(1, nq+1)} - self.e0 = {str(i): 0. for i in range(1, nq+1)} - self.Qp = {'+': {str(i): 0. for i in range(1, nq+1)}, '-':{str(i): 0. for i in range(1, nq+1)}} - - self.basic_forces = np.array([IntForce(self, i) for i in range(1, nq+1)]) - self.basic_deformations = self.basic_forces - - @property - def force_keys(self): - return [self.tag+'_'+key for key in self.rel] - - - def setHistoryParameter ( self, name, val ): - self.current[name] = val - - - def getHistoryParameter ( self, name ): - return self.history[name] - - - def commitHistory ( self ): - self.history = self.current.copy() - self.current = {} - - if hasattr( self , "mat" ): - self.mat.commitHistory() - - def v0_vector(self): - return np.array([0]*self.ndf*self.nn) - - def pw_vector(self): - if all(self.w.values())==0.0: - return np.array([0.0]*self.ndf*self.nn) - -class PolyRod(Element): - nv = 1 - nn = 2 - ndm = 1 - ndf = 1 # number of dofs at each node - force_dict = {'1':0} - - def __init__(self,tag, nodes, E, A): - super().__init__(self.ndf, self.ndm, self.force_dict, nodes) - self.tag = tag - self.E = E - self.A = A - self.q = {'1':0} - self.v = {'1':0} - self.w = {'1':0.0} - - if isinstance(self.E,float): - self.E = Polynomial([self.E]) - - if isinstance(self.A,float): - self.A = Polynomial([self.A]) - - def N(self): - L = self.L - N1 = Polynomial([1,-1/L]) - N2 = Polynomial([0, 1/L]) - return np.array([[N1],[N2]]) - - def B(self): - N = self.N() - return np.array([[Polynomial.deriv(Ni,1) for Ni in row] for row in N]) - - def k_matrix(self): - E = self.E - A = self.A - - L = self.L - B = self.B() - k = Structural_Matrix([ - [quad(E*A*(B@B.T)[i,j],0,L)[0] for j in range(2)] for i in range(2) - ]) - - # Metadata - k.tag = 'k' - k.row_data = k.column_data = ['u_'+str(int(dof)) for dof in self.dofs] - return k - - def ke_matrix(self): - return self.k_matrix() - - def pw_vector(self): - L = self.L - pw = self.w['1'] - N = self.N() - if isinstance(pw,np.polynomial.Polynomial) or isinstance(pw,float): - # p1 = -quad(N[0]*pw,0,L)[0] - # p2 = -quad(N[1]*pw,0,L)[0] - p = np.array([[-quad(Ni*pw,0,L)[0] for Ni in row] for row in N]) - else: - print('Unsupported load vector pw') - - return p - - ## Post Processing - def localize(self,U_vector): - dofs = [int(dof) for dof in self.dofs] - return np.array([U_vector.get(str(dof)) for dof in dofs]) - - def displx(self,U_vector): - u = self.localize(U_vector) - N = self.N() - return N[0,0]*u[0] + N[1,0]*u[1] - - def strainx(self,U_vector): - dofs = [int(dof) for dof in self.dofs] - u = np.array([U_vector.get(str(dof)) for dof in dofs]) - B = self.B() - return B[0,0]*u[0] + B[1,0]*u[1] - - def iforcex(self,U_vector): - """Resisting force vector""" - dofs = [int(dof) for dof in self.dofs] - u = np.array([U_vector.get(str(dof)) for dof in dofs]) - B = self.B() - P = self.E*self.A*(B[0,0]*u[0] + B[1,0]*u[1]) - return P - -class Truss(Element): - nv = 1 - nn = 2 - nq = 1 - ndm = 2 - ndf = 2 # number of dofs at each node - force_dict = {'1':0} - Qpl = np.zeros((2,nq)) - - def __init__(self, tag, iNode, jNode, E=None, A=None, geo='lin',properties=None): - super().__init__(self.ndf, self.ndm, self.force_dict, [iNode, jNode]) - if isinstance(properties,dict): - E, A = properties['E'], properties['A'] - self.type = '2D truss' - self.tag = tag - self.E = E - self.A = A - self.geo=geo - self.Q = np.array([0.0]) - - self.q = {'1':0} - self.v = {'1':0} - - self.w = {'1':0.0} - - def __repr__(self): - return 'truss-{}'.format(self.tag) - - def N(self): - L = self.L - N1 = Polynomial([1,-1/L]) - N2 = Polynomial([0, 1/L]) - return np.array([N1,N2]) - - def B(self): - L = self.L - B1 = self.N()[0].deriv(1) - B2 = self.N()[1].deriv(1) - return np.array([B1,B2]) - - def v0_vector(self): - EA = self.E*self.A - L = self.L - e0 = self.e0 - q0 = self.q0 - w = self.w - v0 = np.array([e0['1']*L]) - v0 += [q0['1']*L/EA] - v0 += [w['1']*L*L/(2*EA)] - return v0 - - def q0_vector(self): - EA = self.E*self.A - e0 = self.e0['1'] - q0 = self.q0['1'] - q0 = q0 - EA*e0 - return [q0] - - def pw_vector(self): - L = self.L - pw = self.w['1'] - N = self.N() - if isinstance(pw,np.polynomial.Polynomial) or isinstance(pw,float): - p1 = -quad(N[0]*pw,0,L)[0] - p2 = -quad(N[1]*pw,0,L)[0] - else: - print('Unsupported load vector pw') - - return np.array([p1,0,p2,0]) - - def ke_matrix(self): - ag = self.ag() - k = self.k_matrix() - return ag.T@(k*ag) - - def kg_matrix(self, N): - """return element local stiffness Matrix""" - E = self.E - A = self.A - L = self.L - k = Structural_Matrix([N/L]) - # Metadata - k.tag = 'kg' - k.row_data = ['q_'+ key for key in self.q0.keys()] - k.column_data = ['v_' + key for key in self.e0.keys()] - k.c_cidx = k.c_ridx = [int(key)-1 for key in self.rel.keys() if self.rel[key]] - return k - - def k_matrix(self): - """return element local stiffness Matrix""" - - E = self.E - A = self.A - L = self.L - k = Structural_Matrix([E*A/L]) - # Metadata - k.tag = 'k' - k.row_data = ['q_'+ key for key in self.q0.keys()] - k.column_data = ['v_' + key for key in self.e0.keys()] - k.c_cidx = k.c_ridx = [int(key)-1 for key in self.rel.keys() if self.rel[key]] - return k - - def bg_matrix(self): - """return element static matrix, :math:`\\mathbf{b}_g`""" - - cs = self.cs - sn = self.sn - bg = np.array([[-cs], - [-sn], - [ cs], - [ sn]]) - return bg - - def f_matrix(self, Roption=True): - """return element flexibility matrix, :math:`\\mathbf{f}`""" - - A = self.A - E = self.E - L = self.L - f = Structural_Matrix([L/(E*A)]) - - if Roption: - pass - - # Metadata - f.tag = 'f' - f.column_data = ['q_' + key for key in self.q0.keys()] - f.row_data = ['v_'+ key for key in self.e0.keys()] - f.c_cidx = f.c_ridx = [int(key)-1 for key in self.rel.keys() if self.rel[key]] - - return f - - def ag(self): - cs = self.cs - sn = self.sn - ag = np.array([[-cs,-sn , cs, sn],]) - return ag - - - def GLstrain(self): - Li = self.Li - L0 = self.L0 - E_GL = (Li**2 - L0**2) / (2*L0**2) - return E_GL - - def iGLstrain(self): - """incremental Green-Lagrange strain""" - L = self.L - Li = self.Li - E_GL = (L**2 - Li**2) / (2*Li**2) - return E_GL - -class TaperedTruss(Truss): - def k_matrix(self): - if isinstance(self.E,float): - E = Polynomial([self.E]) - else: - E = self.E - if isinstance(self.A,float): - A = Polynomial([self.A]) - else: - A = self.A - A = self.A - L = self.L - B = self.B() - # ke = np.zeros((self.ndf,self.ndf)) - # for i in range(self.ndf): - # for j in range(self.ndf): - # f = lambda x: B[i](x)*E(x)*A(x)*B[j](x) - # ke[i,j] = quad(f,0,L)[0] - - f = lambda x: B[0](x)*E(x)*A(x)*B[0](x) - k = Structural_Matrix([quad(f,0,L)[0]]) - - # Metadata - k.tag = 'k' - k.row_data = ['q_'+ key for key in self.q0.keys()] - k.column_data = ['v_' + key for key in self.e0.keys()] - k.c_cidx = k.c_ridx = [int(key)-1 for key in self.rel.keys() if self.rel[key]] - return k - - def q0_vector(self): - # EA = self.E*self.A() - # e0 = self.e0['1'] - # q0 = self.q0['1'] - # q0 = q0 - EA*e0 - return [0] - - def v0_vector(self): - return [0] - -class Truss3D(Element): - ndf = 3 - ndm = 3 - force_dict = {'1':0} - def __init__(self, tag, iNode, jNode, A, E): - super().__init__(self.ndf, self.ndm, self.force_dict, [iNode,jNode]) - self.type = '2D truss' - self.tag = tag - self.A = A - self.E = E - - self.q = {'1':0} - self.v = {'1':0} - - def __repr__(self): - return 'tr-{}'.format(self.tag) - - def bg_matrix(self): - """return element static matrix, bg - pp. 57""" - cs = self.cs - sn = self.sn - cz = self.cz - bg = np.array([[-cs], - [-sn], - [-cz], - [ cs], - [ sn], - [ cz]]) - return bg - - -class Frame(Element): - """linear 2D Euler-Bernouli frame element""" - - nv = 3 - nn = 2 - nq = 3 - ndm = 2 - ndf = 3 - force_dict = {'1':0, '2':0, '3': 0} - Qpl = np.zeros((2,nq)) - - def __init__(self, tag, iNode, jNode, E=None, A=None, I=None,properties=None): - super().__init__(self.ndf, self.ndm, self.force_dict, [iNode,jNode]) - if isinstance(properties,dict): - E, A, I = properties['E'], properties['A'], properties['I'] - self.type = '2D beam' - self.tag = tag - self.E = E - self.A = A - self.I = I - self.q = {'1':0, '2':0, '3': 0} # basic element force - self.v = {'1':0, '2':0, '3': 0} - self.w = {'x':0, 'y':0} #': "uniform element loads in x and y", - - def __repr__(self): - return 'el-{}'.format(self.tag) - - def elastic_curve(self, x, end_rotations, scale=10, global_coord=False): - n = len(x) - L = self.L - N1 = 1-3*(x/L)**2+2*(x/L)**3 - N2 = L*(x/L-2*(x/L)**2+(x/L)**3) - N3 = 3*(x/L)**2-2*(x/L)**3 - N4 = L*((x/L)**3-(x/L)**2) - vi = end_rotations[0] - vj = end_rotations[1] - y = np.array(vi*N2+vj*N4)*scale - xy = np.concatenate(([x],[y])) - if global_coord: - x0 = self.nodes[0].x - y0 = self.nodes[0].y - Rz = self.Rz_matrix()[0:2,0:2] - xy = Rz@xy + [[x0]*n,[y0]*n] - return xy - - @property - def enq(self): - """element number of forces, considering deprecation""" - # output like [1, 1, 1] - return [1 for x in self.q] - - def ag(self): - cs = self.cs - sn = self.sn - L = self.L - ag = np.array([[ -cs , -sn , 0, cs , sn , 0], - [-sn/L, cs/L, 1, sn/L, -cs/L, 0], - [-sn/L, cs/L, 0, sn/L, -cs/L, 1]]) - - if self.dofs[0] == self.dofs[3] or self.dofs[1] == self.dofs[4]: - ag[0,:] = [0.0]*ag.shape[1] - return ag - - def ah(self): - MR = [1 if x else 0 for x in self.rel.values()] - ah = np.array([[1-MR[0], 0 , 0 ], - [ 0 , 1-MR[1] , -0.5*(1-MR[2])*MR[1]], - [ 0 , -0.5*(1-MR[1])*MR[2], 1-MR[2] ]]) - return ah - - def k_matrix(self): - """return element local stiffness Matrix""" - E = self.E - A = self.A - I = self.I - L = self.L - EI = E*I - ah = self.ah() - k = np.array([[E*A/L, 0 , 0 ], - [ 0 , 4*EI/L , 2*EI/L], - [ 0 , 2*EI/L , 4*EI/L]]) - k = ah.T @ k @ ah - - # Assemble matrix metadata - k = Structural_Matrix(k) - k.tag = 'k' - k.row_data = ['q_'+ key for key in self.q0.keys()] - k.column_data = ['v_' + key for key in self.v0.keys()] - k.c_cidx = k.c_ridx = [int(key)-1 for key in self.rel.keys() if self.rel[key]] - return k - - def f_matrix(self, Roption=False): - """Flexibility matrix of an element. - - """ - EA = self.E*self.A - EI = self.E*self.I - L = self.L - f = Structural_Matrix([ - [L/EA, 0 , 0 ], - [ 0 , L/(3*EI), -L/(6*EI)], - [ 0 , -L/(6*EI), L/(3*EI)]]) - ide = set(range(3)) - - if Roption: - if self.rel["2"]: - f[:,1] = [0.0]*f.shape[0] - f[1,:] = [0.0]*f.shape[1] - - if self.rel["3"]: - f[:,2] = [0.0]*f.shape[0] - f[2,:] = [0.0]*f.shape[1] - - # Define matrix metadata - f.tag = 'f' - f.column_data = ['q_'+ key for key in self.q0.keys()] - f.row_data = ['v_' + key for key in self.v0.keys()] - f.c_cidx = f.c_ridx = [int(key)-1 for key in self.rel.keys() if self.rel[key]] - return f - - def ke_matrix(self): - """return element global stiffness Matrix""" - - k = self.k_matrix() - ag = self.ag() - - ke = ag.T @ k @ ag - ke = Structural_Matrix(ke) - ke.row_data = ke.column_data = ['u_'+str(int(dof)) for dof in self.dofs] - return ke - - def bg_matrix(self, Roption=False): - """return element global static matrix, bg""" - - cs = self.cs - sn = self.sn - L = self.L - # x ri rj Global - bg = np.array([[-cs, -sn/L, -sn/L], # x - [-sn, cs/L, cs/L], # y - [0.0, 1.0, 0.0], # rz - [ cs, sn/L, sn/L], - [ sn, -cs/L, -cs/L], - [0.0, 0.0, 1.0]]) - if Roption: - if self.rel['2']: - bg[:,1] = [0.0]*bg.shape[0] - - if self.rel['3']: - bg[:,2] = [0.0]*bg.shape[0] - - if self.dofs[0] == self.dofs[3] or self.dofs[1] == self.dofs[4]: - bg[:,0] = [0.0]*bg.shape[0] - # bg = np.delete(bg, 0, axis=1) - - bg = Structural_Matrix(bg) - bg.tag = 'b' - bg.column_data = ['x', 'ri', 'rj'] - bg.row_data = ['x', 'y', 'rz'] - bg.c_cidx = bg.c_ridx = [int(key)-1 for key in self.rel.keys() if self.rel[key]] - - - return bg - - def v0_vector(self): - EA = self.E*self.A - EI = self.E*self.I - L = self.L - e0 = self.e0 - w = self.w - v0 = np.array([e0['1']*L, -e0['2']*L/2, e0['2']*L/2]) - v0 += np.array([w['x']*L*L/(2*EA), w['y']*L**3/(24*EI), -w['y']*L**3/(24*EI)]) - v0 = Structural_Matrix(v0) - v0.tag = 'v' - v0.column_data = ['v_0'] - v0.row_data = ['axial', 'i-rotation', 'j-rotation'] - v0.c_cidx = False - v0.c_ridx = [int(key)-1 for key in self.rel.keys() if self.rel[key]] - return v0 - - def q0_vector(self): - L = self.L - E = self.E - A = self.A - I = self.I - e0= self.e0 - w = self.w - q0 = np.array([-e0['1']*E*A, +e0['2']*E*I, -e0['2']*E*I]) - if self.rel['2']: - q0[1] *= 0 - q0[1] *= 1.5 - if self.rel['3']: - q0[2] *= 0 - q0[2] *= 1.5 - - q0[1] += -w['y']*L**2/12*((not self.rel['2']) and (not self.rel['3'])) - w['y']*L**2/8*(self.rel['3']) - q0[2] += +w['y']*L**2/12*((not self.rel['2']) and (not self.rel['3'])) + w['y']*L**2/8*(self.rel['2']) - - # Metadata - q0 = Structural_Matrix(q0) - q0.tag = 'q' - q0.column_data = ['q_0'] - q0.row_data = ['axial', 'M_i', 'M_j'] - q0.c_cidx = False - q0.c_ridx = [int(key)-1 for key in self.rel.keys() if self.rel[key]] - - return q0 - - def κ(self, k, state=None, form=None): - """Defines element curvature and calculates the resulting end deformations. - - """ - if form ==None: form = 'uniform' - if state==None: state = -1 - L = self.L - table = { - 'uniform': [-0.5*k*L, 0.5*k*L], - 'ilinear': [-k/3*L, k/6*L], - 'jlinear': [-k/6*L, k/2*L], - 'parabolic': [-k/3*L, k/3*L], - } - self.nodes[0].model.states[state]['v0'][self.tag]['2'] = table[form][0] - self.nodes[0].model.states[state]['v0'][self.tag]["3"]= table[form][1] - return - - -class PrismFrame(Frame): - axial_force = 0.0 - - def k_matrix_exact(self, axial_force=0.0): - P = axial_force - - L = self.L - EI = self.E*self.I - k = np.sqrt(abs(P/EI)) - kL = k*L - d = 2*(1-np.cos(kL))-(kL*np.sin(kL)) - ka = self.E*self.A / L - kfv = (EI/(L**3))*(kL**3)*np.sin(kL)/d - kmv = (EI/(L**2))*((kL**2)*(1-np.cos(kL))/d) - kft = kmv - kmt =(EI/L)*((kL*(np.sin(kL)-kL*np.cos(kL)))/d) - kmth = (EI/L)*((kL*(kL-np.sin(kL)))/d) - - eltk = np.array([ - [ ka, 0.0, 0.0, -ka, 0.0, 0.0], - [ 0.0, kfv, kft, 0.0, -kfv, kft], - [ 0.0, kmv, kmt, 0.0, -kmv, kmth], - [ -ka, 0.0, 0.0, ka, 0.0, 0.0], - [ 0.0, -kfv, -kft, 0.0, kfv, -kft], - [ 0.0, kft, kmth, 0.0, -kft, kmt]]) - - # Assemble matrix metadata - k = Structural_Matrix(eltk) - k.tag = 'k' - k.row_data = ['q_'+ key for key in self.q0.keys()]*2 - k.column_data = ['v_' + key for key in self.v0.keys()]*2 - k.c_cidx = k.c_ridx = [int(key)-1 for key in self.rel.keys() if self.rel[key]] - return k - - def k_matrix(self, axial_force=None): - """Evaluates the 2D beam-column stiffness matrix using a - power series expansion - """ - if axial_force is None: P = self.axial_force - else: P = axial_force - L = self.L - EI = self.E*self.I - k = np.sqrt(P/EI) - kL = k*L - ka = self.E*self.A / L - kfv = 12*EI/(L**3)-6*EI/(5*L)*k**2-EI*L/700*k**4 - kmv = 6*EI/L**2 - EI/10*k**2 - EI*L**2/1400*k**4 - kft = kmv - kmt = 4*EI/L-2*EI*L/15*k**2 - 11*EI/6300*L**3*k**4 - kmth = 2*EI/L+EI*L/30*k**2 + 13*EI*L**3*k**4/12600 - - eltk = np.array([ - [ ka, 0.0, -0.0, -ka, 0.0, 0.0], - [ 0.0, kfv, kft, -0.0, -kfv, kft], - [ 0.0, kmv, kmt, -0.0, -kmv, kmth], - [ -ka, 0.0, -0.0, ka, 0.0, 0.0], - [ 0.0, -kfv, -kft, -0.0, kfv, -kft], - [ 0.0, kft, kmth, -0.0, -kft, kmt]]) - - # hinges - # cols = rows = [int(key)-1 for key in self.rel.keys() if self.rel[key]] - # f = np.linalg.inv(eltk) - - # for row in rows: - # f[row, :] = [0.0]*6 - # for col in cols: - # f[:, col] = [0.0]*6 - # k = np.linalg.inv(f) - - - # Matrix metadata - k = Structural_Matrix(eltk) - k.tag = 'k' - k.row_data = ['q_'+ key for key in self.q0.keys()]*2 - k.column_data = ['v_' + key for key in self.v0.keys()]*2 - k.c_cidx = k.c_ridx = [int(key)-1 for key in self.rel.keys() if self.rel[key]] - return k - - def ke_matrix(self, axial_force=None): - T = self.Rz_matrix() - k = self.k_matrix(axial_force) - ke = T.T @ k @ T - return ke - -class Beam3d(Element): - nv = 6 - nn = 2 - ndm = 3 - ndf = 6 - force_dict = {'1':0, '2':0, '3': 0} - - def __init__(self, tag, iNode, jNode, E, A, Iy, Iz, Gs, Kv): - super().__init__(self.ndf, self.ndm, self.force_dict, [iNode,jNode]) - self.type = '3D beam' - self.tag = tag - self.E = E - self.A = A - self.Gs = Gs - self.Iy = Iy - self.Iz = Iz - self.Kv = Kv - # self.q = {'1':0, '2':0, '3': 0} # basic element force - # self.v = {'1':0, '2':0, '3': 0} - - self.w = {'x':0, 'y':0, 'z': 0} #': "uniform element loads in x and y", - - - def __repr__(self): - return 'el-{}'.format(self.tag) - - def k_matrix(self): - """ - Calculate the stiffness matrix for a 3D elastic Bernoulli - beam element. - - Parameters - -------------------------------------------- - E: Young's modulus - G: Shear modulus - A: Cross section area - Iy: Moment of inertia, local y-axis - Iz: Moment of inertia, local z-axis - Kv: Saint-Venant's torsion constant - - Returns: - Kle local beam stiffness matrix (12 x 12) - """ - - L,E,Gs,A,Iy,Iz,Kv = self.L, self.E, self.Gs, self.A, self.Iy, self.Iz, self.Kv - - a = E*A/L - b = 12*E*Iz/L**3 - c = 6*E*Iz/L**2 - d = 12*E*Iy/L**3 - e = 6*E*Iy/L**2 - f = Gs*Kv/L - g = 2*E*Iy/L - h = 2*E*Iz/L - - Kle = np.mat([ - [ a, 0, 0, 0, 0, 0, -a, 0, 0, 0, 0, 0 ], - [ 0, b, 0, 0, 0, c, 0, -b, 0, 0, 0, c ], - [ 0, 0, d, 0, -e, 0, 0, 0, -d, 0, -e, 0 ], - [ 0, 0, 0, f, 0, 0, 0, 0, 0, -f, 0, 0 ], - [ 0, 0, -e, 0, 2*g, 0, 0, 0, e, 0, g, 0 ], - [ 0, c, 0, 0, 0, 2*h, 0, -c, 0, 0, 0, h ], - [-a, 0, 0, 0, 0, 0, a, 0, 0, 0, 0, 0 ], - [ 0, -b, 0, 0, 0, -c, 0, b, 0, 0, 0, -c ], - [ 0, 0, -d, 0, e, 0, 0, 0, d, 0, e, 0 ], - [ 0, 0, 0, -f, 0, 0, 0, 0, 0, f, 0, 0 ], - [ 0, 0, -e, 0, g, 0, 0, 0, e, 0, 2*g, 0 ], - [ 0, c, 0, 0, 0, h, 0, -c, 0, 0, 0, 2*h ]]) - - # fle = L/2*np.mat([qx, qy, qz, qw, -qz*L/6, qy*L/6, qx, qy, qz, qw, qz*L/6, -qy*L/6]).T - - # n2 = np.array([0.,0.,0.]) - # n2[0] = n3[1]*n1[2]-n3[2]*n1[1] - # n2[1] = -n1[2]*n3[0]+n1[0]*n3[2] - # n2[2] = n3[0]*n1[1]-n1[0]*n3[1] - - #An = np.append([n1,n2],[n3],0) - - # G = np.mat([ - # [ n1[0], n1[1], n1[2], 0, 0, 0, 0, 0, 0, 0, 0, 0 ], - # [ n2[0], n2[1], n2[2], 0, 0, 0, 0, 0, 0, 0, 0, 0 ], - # [ n3[0], n3[1], n3[2], 0, 0, 0, 0, 0, 0, 0, 0, 0 ], - # [ 0, 0, 0, n1[0], n1[1], n1[2], 0, 0, 0, 0, 0, 0 ], - # [ 0, 0, 0, n2[0], n2[1], n2[2], 0, 0, 0, 0, 0, 0 ], - # [ 0, 0, 0, n3[0], n3[1], n3[2], 0, 0, 0, 0, 0, 0 ], - # [ 0, 0, 0, 0, 0, 0, n1[0], n1[1], n1[2], 0, 0, 0 ], - # [ 0, 0, 0, 0, 0, 0, n2[0], n2[1], n2[2], 0, 0, 0 ], - # [ 0, 0, 0, 0, 0, 0, n3[0], n3[1], n3[2], 0, 0, 0 ], - # [ 0, 0, 0, 0, 0, 0, 0, 0, 0, n1[0], n1[1], n1[2]], - # [ 0, 0, 0, 0, 0, 0, 0, 0, 0, n2[0], n2[1], n2[2]], - # [ 0, 0, 0, 0, 0, 0, 0, 0, 0, n3[0], n3[1], n3[2]] - # ]) - - # Ke = G.T*Kle*G - # fe = G.T*fle - return Kle - # if eq == None: - # return Kle - # else: - # return Kle,fe - - - def ke_matrix(self): - Kle = self.k_matrix() - ex = [self.nodes[0].x, self.nodes[1].x] - ey = [self.nodes[0].y, self.nodes[1].y] - ez = [self.nodes[0].z, self.nodes[1].z] - b = np.mat([ - [ex[1]-ex[0]], - [ey[1]-ey[0]], - [ez[1]-ez[0]]]) - - L = self.L - n1 = np.asarray(b.T/L).reshape(3,) - - eo = np.asmatrix(eo) - lc = np.asscalar(np.sqrt(eo*eo.T)) - n3 = np.asarray(eo/lc).reshape(3,) - - - n2 = np.array([0.,0.,0.]) - n2[0] = n3[1]*n1[2]-n3[2]*n1[1] - n2[1] = -n1[2]*n3[0]+n1[0]*n3[2] - n2[2] = n3[0]*n1[1]-n1[0]*n3[1] - - T = np.mat([ - [ n1[0], n1[1], n1[2], 0, 0, 0, 0, 0, 0, 0, 0, 0 ], - [ n2[0], n2[1], n2[2], 0, 0, 0, 0, 0, 0, 0, 0, 0 ], - [ n3[0], n3[1], n3[2], 0, 0, 0, 0, 0, 0, 0, 0, 0 ], - [ 0, 0, 0, n1[0], n1[1], n1[2], 0, 0, 0, 0, 0, 0 ], - [ 0, 0, 0, n2[0], n2[1], n2[2], 0, 0, 0, 0, 0, 0 ], - [ 0, 0, 0, n3[0], n3[1], n3[2], 0, 0, 0, 0, 0, 0 ], - [ 0, 0, 0, 0, 0, 0, n1[0], n1[1], n1[2], 0, 0, 0 ], - [ 0, 0, 0, 0, 0, 0, n2[0], n2[1], n2[2], 0, 0, 0 ], - [ 0, 0, 0, 0, 0, 0, n3[0], n3[1], n3[2], 0, 0, 0 ], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, n1[0], n1[1], n1[2]], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, n2[0], n2[1], n2[2]], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, n3[0], n3[1], n3[2]] - ]) - Ke = T.T*Kle*T - return Ke - - diff --git a/ops/opensees/emit/__init__.py b/SRC/opensees/emit/__init__.py similarity index 100% rename from ops/opensees/emit/__init__.py rename to SRC/opensees/emit/__init__.py diff --git a/ops/opensees/emit/__main__.py b/SRC/opensees/emit/__main__.py similarity index 100% rename from ops/opensees/emit/__main__.py rename to SRC/opensees/emit/__main__.py diff --git a/ops/opensees/emit/apidoc.py b/SRC/opensees/emit/apidoc.py similarity index 100% rename from ops/opensees/emit/apidoc.py rename to SRC/opensees/emit/apidoc.py diff --git a/ops/opensees/emit/emit_json.py b/SRC/opensees/emit/emit_json.py similarity index 100% rename from ops/opensees/emit/emit_json.py rename to SRC/opensees/emit/emit_json.py diff --git a/ops/opensees/emit/emitter.py b/SRC/opensees/emit/emitter.py similarity index 100% rename from ops/opensees/emit/emitter.py rename to SRC/opensees/emit/emitter.py diff --git a/ops/opensees/emit/fedeas.py b/SRC/opensees/emit/fedeas.py similarity index 100% rename from ops/opensees/emit/fedeas.py rename to SRC/opensees/emit/fedeas.py diff --git a/ops/opensees/emit/mesh.py b/SRC/opensees/emit/mesh.py similarity index 100% rename from ops/opensees/emit/mesh.py rename to SRC/opensees/emit/mesh.py diff --git a/ops/opensees/emit/opensees.py b/SRC/opensees/emit/opensees.py similarity index 100% rename from ops/opensees/emit/opensees.py rename to SRC/opensees/emit/opensees.py diff --git a/ops/opensees/emit/registry.py b/SRC/opensees/emit/registry.py similarity index 100% rename from ops/opensees/emit/registry.py rename to SRC/opensees/emit/registry.py diff --git a/ops/opensees/emit/writer.py b/SRC/opensees/emit/writer.py similarity index 100% rename from ops/opensees/emit/writer.py rename to SRC/opensees/emit/writer.py diff --git a/ops/opensees/library/__init__.py b/SRC/opensees/library/__init__.py similarity index 100% rename from ops/opensees/library/__init__.py rename to SRC/opensees/library/__init__.py diff --git a/ops/opensees/library/ast.py b/SRC/opensees/library/ast.py similarity index 100% rename from ops/opensees/library/ast.py rename to SRC/opensees/library/ast.py diff --git a/ops/opensees/library/material.py b/SRC/opensees/library/material.py similarity index 100% rename from ops/opensees/library/material.py rename to SRC/opensees/library/material.py diff --git a/ops/opensees/library/model.py b/SRC/opensees/library/model.py similarity index 100% rename from ops/opensees/library/model.py rename to SRC/opensees/library/model.py diff --git a/ops/opensees/library/obj.py b/SRC/opensees/library/obj.py similarity index 100% rename from ops/opensees/library/obj.py rename to SRC/opensees/library/obj.py diff --git a/ops/opensees/library/pattern.py b/SRC/opensees/library/pattern.py similarity index 100% rename from ops/opensees/library/pattern.py rename to SRC/opensees/library/pattern.py diff --git a/SRC/opensees/material.py b/SRC/opensees/material.py deleted file mode 100644 index 8b13789179..0000000000 --- a/SRC/opensees/material.py +++ /dev/null @@ -1 +0,0 @@ - diff --git a/SRC/opensees/matrices.py b/SRC/opensees/matrices.py deleted file mode 100644 index dc52142618..0000000000 --- a/SRC/opensees/matrices.py +++ /dev/null @@ -1,2333 +0,0 @@ -""" -Matrices. (:mod:`matrices`) -===================================================== - -.. currentmodule:: matrices - -This module provides functions and classes for constructing various -structural analysis matrices. Matrices are generally constructed with -functions that are defined outside of classes (as opposed to using the -__init__ method of a class) in an effort to present basic structural analysis -concepts with minimal clutter from computer science concepts which may -be unfamiliar to engineers. - -""" - -import numpy as np -import pandas as pd -import scipy.linalg -import matplotlib.pyplot as plt - - -import numpy as np -import matplotlib.pyplot as plt -import copy - - -def _create_model(model): - pass - -def static_matrix(model): - pass - -def strain_matrix(model): - pass - -##<****************************Model objects************************************ -# Model objects/classes -# These should be created using the methods above -#****************************************************************************** -class _Node: - def __init__(self, model, tag: str, ndf, xyz): - - self.xyz = np.array([xi for xi in xyz if xi is not None]) - - self.tag = tag - self.xyz0 = self.xyz # coordinates in base configuration (unstrained, not necessarily unstressed). - self.xyzi = self.xyz # coordinates in reference configuration. - - self.x0: float = xyz[0] # x-coordinate in base configuration (unstrained, not necessarily unstressed). - self.y0: float = xyz[1] # y-coordinate in base configuration (unstrained, not necessarily unstressed). - self.z0: float = xyz[2] # z-coordinate in base configuration (unstrained, not necessarily unstressed). - - - self.x: float = xyz[0] - self.y: float = xyz[1] - self.z: float = xyz[2] - - - self.rxns = [0]*ndf - self.model = model - self.elems = [] - - self.p = {dof:0.0 for dof in model.ddof} - - def __repr__(self): - return 'nd-{}'.format(self.tag) - - def p_vector(self): - return np.array(list(self.p.values())) - - @property - def dofs(self): - # if self.model.DOF == None: self.model.numDOF() - idx = self.model.nodes.index(self) - return np.array(self.model.DOF[idx],dtype=int) - - -class Rxn: - def __init__(self, node, dirn): - self.node = node - self.dirn = dirn - - def __repr__(self): - return 'rxn-{}'.format(self.dirn) - - -class Hinge: - def __init__(self, elem, node): - self.elem = elem - self.node = node - - -class Model: - def __init__(self, ndm, ndf): - """Basic model object - - Parameters - ----------- - ndm: - number of model dimensions - ndf: - number of degrees of freedom (dofs) at each node - - """ - self.ndf: int = ndf - self.ndm: int = ndm - self.DOF: list = None - - # Define DOF list indexing - if ndm == 1: - self.prob_type = '1d' - self.ddof: dict = {'x': 0} # Degrees of freedom at each node - - if ndf == 2: - self.prob_type = '2d-truss' - self.ddof: dict = { 'x': 0, 'y': 1} # Degrees of freedom - elif ndm == 2 and ndf ==3: - self.prob_type = '2d-frame' - self.ddof: dict = { 'x': 0, 'y': 1, 'rz':2} - elif ndm == 3 and ndf ==3: - self.prob_type = '3d-truss' - self.ddof: dict = { 'x': 0, 'y': 1, 'z':2} - elif ndm == 3 and ndf ==6: - self.prob_type = '3d-frame' - self.ddof: dict = { 'x': 0, 'y': 1, 'z':2, 'rx':3, 'ry':4, 'rz':5} - - # model inventory lists - self.elems: list = [] - self.nodes: list = [] - self.rxns: list = [] - self.hinges: list = [] - self.iforces: list = [] - self.loads: list = [] - self.states: list = [] - self.redundants: list = [] - - # model inventory dictionaries - self.delems: dict = {} - self.dnodes: dict = {} - self.dxsecs: dict = {} - self.dhinges: dict = {} - self.dstates: dict = {} - self.dxsecs: dict = {} - self.materials: dict = {} - self.xsecs: dict = {} - self.dredundants: dict = {} - - # Initialize default material/section properties - self.material('default', 1.0) - self.xsection('default', 1.0, 1.0) - - @property - def rel(self): - return [rel for elem in self.elems for rel in elem.rel.values()] - - @property - def nn(self) -> int: - """return number of nodes in model""" - return len(self.nodes) - - @property - def nr(self) -> int: - """return number of constrained dofs in model""" - return len(self.rxns) - - @property - def ne(self) -> int: - """return number of elements in model""" - return len(self.elems) - - @property - def nQ(self): - f = 0 - for elem in self.elems: - f += len(elem.basic_forces) - return f - - @property - def nq(self): - f = [] - for elem in self.elems: - f.append(sum([1 for x in elem.q])) - return f - - @property - def nv(self): - lst = [] - for elem in self.elems: - lst.append(sum([1 for x in elem.v])) - return lst - - @property - def nf(self) -> int: - x = self.nt - self.nr - return x - - @property - def nt(self) -> int: - return self.ndf*self.nn - - @property - def rdofs(self): - """Return list of restrained dofs""" - DOF = self.DOF - - return [] - - @property - def NOS(self) -> int: - nf = self.nf - nq = sum(self.nq) - return nq - nf - - @property - def basic_forces(self): - return np.array([q for elem in self.elems for q in elem.basic_forces ]) - - @property - def rdnt_forces(self): - cforces = self.cforces - return np.array([q for q in cforces if q.redundant ]) - - @property - def cforces(self): - return np.array([q for elem in self.elems for q in elem.basic_forces if not q.rel]) - - @property - def eforces(self): - """Return array of elastic element forces""" - return np.array([q for elem in self.elems for q in elem.basic_forces if (q.plastic_event is None)]) - - @property - def idx_c(self): - cforces = self.cforces - forces = self.basic_forces - idx_c = np.where(np.isin(forces,cforces))[0] - return idx_c - - @property - def idx_e(self): - """return indices of elastic basic (not plastic) forces""" - cforces = self.cforces - eforces = self.eforces - idx_e = np.where(np.isin(cforces,eforces))[0] - return idx_e - - @property - def idx_f(self): - return np.arange(0,self.nf) - @property - def idxx_f(self): - return None - - @property - def idx_i(self): - rdts = self.rdnt_forces - forces = self.basic_forces - idx_i = np.where(np.logical_not(np.isin(forces, rdts)))[0] - return idx_i - - @property - def idx_x(self): - rdts = self.rdnt_forces - forces = self.basic_forces - idx_x = np.where(np.isin(forces, rdts))[0] - return idx_x - - def node(self, tag: str, x: float, y=None, z=None, mass: float=None): - newNode = _Node(self, tag, self.ndf, [x, y, z], mass) - self.nodes.append(newNode) - self.dnodes.update({newNode.tag : newNode}) - return newNode - - def state(self, method="Linear"): - if self.DOF==None: self.numDOF() - newState = State(self, method) - - ElemTypes = {type(elem) for elem in self.elems} - StateVars = {key for elem in ElemTypes for key in elem.stateVars.keys() } - - stateDict = { - var : { - elem.tag : copy.deepcopy(elem.stateVars[var]) - for elem in self.elems if var in elem.stateVars.keys() - } for var in StateVars - } - self.states.append(stateDict) - return stateDict - - def numDOF(self): - crxns = self.ndf*len(self.nodes) - len(self.rxns)+1 - df = 1 - temp = [] - for node in self.nodes: - DOFs = [] - for rxn in node.rxns: - if not(rxn): - DOFs.append(df) - df += 1 - else: - DOFs.append(crxns) - crxns += 1 - temp.append(DOFs) - self.DOF = temp - return self.DOF - - - - def fix(self, node, dirn=["x","y","rz"]): # for dirn enter string (e.g. "x", 'y', 'rz') - if isinstance(dirn,list): - rxns = [] - for df in dirn: - newRxn = Rxn(node, df) - self.rxns.append(newRxn) - rxns.append(newRxn) - node.rxns[self.ddof[df]] = 1 - return rxns - else: - newRxn = Rxn(node, dirn) - self.rxns.append(newRxn) - node.rxns[self.ddof[dirn]] = 1 - return newRxn - - # Other - def material(self, tag: str, E: float): - newMat = Material(tag, E) - self.materials[tag]=newMat - return newMat - - def xsection(self, tag: str, A: float, I: float): - newXSect = XSect(tag, A, I) - self.xsecs[tag] = newXSect - return newXSect - - # Elements - def add_element(self, element): - """Add a general element to model - - Parameters - --------- - element : obj - - """ - - self.elems.append(element) - self.delems.update({element.tag:element}) - - for node in element.nodes: - node.elems.append(element) - - return element - - def add_elements(self, elements): - """Add a general element to model - - Parameters - --------- - element : obj - - """ - for element in elements: - self.elems.append(element) - self.delems.update({element.tag:element}) - for node in element.nodes: - node.elems.append(element) - - return elements - - - def beam(self, tag: str, iNode, jNode, mat=None, sec=None, Qpl=None,**kwds): - """Add a 2D linear Euler-Bernouli beam object to model - - Parameters - --------- - tag : str - string used for identifying object - iNode : ema.Node - node object at element i-end - jNode : ema.Node - node object at element j-end - mat : ema.Material - - sec : ema.Section - - - """ - - if mat is None: - E = kwds["E"] if "E" in kwds else self.materials["default"].E - else: - E = mat.E - - if sec is None: - A = kwds["A"] if "A" in kwds else self.xsecs["default"].A - I = kwds["I"] if "I" in kwds else self.xsecs["default"].I - else: - A = sec.A - I = sec.I - - from ema.elements import Beam - - newElem = Beam(tag, iNode, jNode, E, A, I) - self.elems.append(newElem) - self.delems.update({newElem.tag:newElem}) - # self.connect([iNode, jNode], "Beam") # considering deprecation - iNode.elems.append(newElem) - jNode.elems.append(newElem) - - if Qpl is not None: - newElem.Qpl = np.zeros((3,2)) - newElem.Np = [Qpl[0], Qpl[0]] - newElem.Mp = [[Qpl[1], Qpl[1]],[Qpl[2], Qpl[2]]] - for i, key in enumerate(newElem.Qp['+']): - newElem.Qp['+'][key] = newElem.Qp['-'][key] = Qpl[i] # consider depraction of elem.Qp in favor of elem.Qpl - newElem.Qp['+'][key] = newElem.Qp['-'][key] = Qpl[i] - newElem.Qpl[i,:] = Qpl[i] # <- consider shifting to this format for storing plastic capacities - return newElem - - - def truss(self, tag: str, iNode, jNode, mat=None, xsec=None, Qpl=None,A=None,E=None): - from ema.elements import Truss - - if mat is None: mat = self.materials['default'] - if E is None: E = mat.E - # cross section - if xsec is None: xsec = self.xsecs['default'] - if A is None: A = xsec.A - - - newElem = Truss(tag, iNode, jNode, E, A) - self.delems.update({newElem.tag:newElem}) - self.elems.append(newElem) - iNode.elems.append(newElem) - jNode.elems.append(newElem) - - if Qpl is not None: - newElem.Np = [Qpl[0], Qpl[0]] - newElem.Qp['+']['1'] = newElem.Qp['-']['1'] = Qpl[0] - return newElem - - - def taprod(self, tag: str, iNode, jNode, mat=None, xsec=None, Qpl=None,A=None,E=None): - """Construct a tapered rod element with variable E and A values.""" - from ema.elements import TaperedTruss - - if mat is None: - mat = self.materials['default'] - if E is None: - E = mat.E - # cross section - if xsec is None: xsec = self.xsecs['default'] - if A is None: A = xsec.A - - newElem = TaperedTruss(tag, iNode, jNode, E, A) - self.delems.update({newElem.tag:newElem}) - self.elems.append(newElem) - iNode.elems.append(newElem) - jNode.elems.append(newElem) - - if Qpl is not None: - newElem.Np = [Qpl[0], Qpl[0]] - newElem.Qp['+']['1'] = newElem.Qp['-']['1'] = Qpl[0] - return newElem - - def truss3d(self, tag: str, iNode, jNode, mat=None, xsec=None): - """Add an ema.Truss3d object to model - - Parameters - --------- - - """ - from ema.elements import Truss3D - if mat is None: mat = self.materials['default'] - if xsec is None: xsec = self.xsecs['default'] - newElem = Truss3D(tag, iNode, jNode, mat, xsec) - self.elems.append(newElem) - self.delems.update({newElem.tag:newElem}) - iNode.elems.append(newElem) - jNode.elems.append(newElem) - return newElem - - - def hinge(self, elem, node): - "pin a beam end." - newHinge = Hinge(elem, node) - self.hinges.append(newHinge) - if node == elem.nodes[0]: - elem.rel['2'] = True - elem.q.pop('2') - elem.basic_forces[1].rel = True - elif node == elem.nodes[1]: - elem.rel['3'] = True - elem.q.pop('3') - elem.basic_forces[2].rel = True - else: - raise Exception("element {} is not bound to node {}.".format(elem, node)) - - return newHinge - - - def redundant(self, elem, nature): - newq = IntForce(elem, nature) - elem.red[nature] = True - self.redundants.append(newq) - - -class rModel(Model): - def __init__(self, ndm, ndf): - super().__init__(ndm=2, ndf=3) - self.material('default', 1.0) - self.xsection('default', 1.0, 1.0) - - def isortho(self, elem): - if (abs(elem.cs) == 1.0) or (abs(elem.sn) == 1.0): - return True - else: - return False - - def numdofs(self): - current_rxn = 1 - current_dof = 1 - rxn_ixs = [] - DOFs = [[0, 0, 0] for node in self.nodes] - for i, node in enumerate(self.nodes): - # x-dof - dirn = 0 - if not(node.rxns[dirn]): # node is free - if not(DOFs[i][dirn]): # node unassigned - for elem in node.elems: - if abs(elem.cs) == 1.0: # x-dof coupled to far end - if elem.nodes[0] == node: - far_node = self.nodes.index(elem.nodes[1]) - if elem.nodes[1] == node: - far_node = self.nodes.index(elem.nodes[0]) - - if not(DOFs[far_node][dirn]): # Far node dof unassigned - if not(self.nodes[far_node].rxns[dirn]): # Far node is free - DOFs[far_node][dirn] = current_dof - DOFs[i][dirn] = current_dof - current_dof += 1 - else: # Far node is fixed - DOFs[far_node][dirn] = current_rxn - DOFs[i][dirn] = current_rxn - current_rxn += 1 - rxn_ixs.append( (i,dirn) ) - rxn_ixs.append( (far_node,dirn) ) - else: # Far node dof already assigned - DOFs[i][dirn] = DOFs[far_node][dirn] - if self.nodes[far_node].rxns[dirn]: # Far node is fixed - rxn_ixs.append( (i,dirn) ) - - elif all([abs(elem.cs) != 1.0 for elem in node.elems]): # x-dof free/uncoupled - if not DOFs[i][dirn]: - DOFs[i][dirn] = current_dof - current_dof += 1 - - else: # node is fixed - if not DOFs[i][dirn]: # node is unassigned - DOFs[i][dirn] = current_rxn - current_rxn += 1 - rxn_ixs.append( (i,dirn) ) - - # y-dof - dirn = 1 - if not(node.rxns[dirn]): - if not(DOFs[i][dirn]): - for elem in node.elems: - if abs(elem.sn) == 1.0: - # get far node index - if elem.nodes[0] == node: - far_node = self.nodes.index(elem.nodes[1]) - if elem.nodes[1] == node: - far_node = self.nodes.index(elem.nodes[0]) - - if not(DOFs[far_node][dirn]): - if not(self.nodes[far_node].rxns[dirn]): - DOFs[far_node][dirn] = current_dof - DOFs[i][dirn] = current_dof - current_dof += 1 - else: - DOFs[far_node][dirn] = current_rxn - DOFs[i][dirn] = current_rxn - current_rxn += 1 - rxn_ixs.append( (i,dirn) ) - rxn_ixs.append( (far_node,dirn) ) - else: - DOFs[i][dirn] = DOFs[far_node][dirn] - if self.nodes[far_node].rxns[dirn]: - rxn_ixs.append( (i,dirn) ) - elif all([abs(elem.sn) != 1.0 for elem in node.elems]): - if not(DOFs[i][dirn]): - DOFs[i][dirn] = current_dof - current_dof += 1 - else: - if not(DOFs[i][dirn]): - DOFs[i][dirn] = current_rxn - current_rxn += 1 - rxn_ixs.append( (i,dirn) ) - - # rz-dof - dirn = 2 - if not(node.rxns[2]): - DOFs[i][dirn] = current_dof - current_dof += 1 - else: - DOFs[i][dirn] = current_rxn - current_rxn += 1 - rxn_ixs.append( (i,dirn) ) - - for ids in rxn_ixs: - DOFs[ids[0]][ids[1]] += current_dof - 1 - - return DOFs - - - def numDOF(self): - crxns = self.ndf*len(self.nodes) - len(self.rxns)+1 - df = 1 - temp = [] - for node in self.nodes: - DOFs = [] - for rxn in node.rxns: - if not(rxn): - DOFs.append(df) - df += 1 - else: - DOFs.append(crxns) - crxns += 1 - temp.append(DOFs) - self.DOF = temp - return self.DOF - - @property - def triv_forces(self): - """list of trivial axial forces""" - lst = [] - for elem in self.elems: - if len(elem.basic_forces) > 1: - if elem.dofs[0]==elem.dofs[3] or elem.dofs[1]==elem.dofs[4]: - lst.append(elem.basic_forces[0]) - return np.array(lst) - - # @property - # def basic_forces(self): - # # bmax = self.TrAx_forces - # forces = np.array([q for elem in self.elems for q in elem.basic_forces if not q.rel]) - # return forces - - @property - def cforces(self): - triv = self.triv_forces - arry = np.array([q for elem in self.elems for q in elem.basic_forces if ( - (q.plastic_event is None) and ( - not q in triv) and ( - not q.rel))]) - return arry - - @property - def nr(self): - return len(self.rxns) - - # @property - # def nq(self): - # f = [] - # for elem in self.elems: - # f.append(sum([1 for q in elem.basic_forces if not q.rel and (not q in self.triv_forces)])) - # return f - - @property - def nv(self): - """Returns number of element deformations in model""" - lst = [] - for elem in self.elems: - lst.append(sum([1 for x in elem.v])) - return lst - - - - @property - def fdof(self): - """Return list of free dofs""" - pass - - @property - def nt(self): - nt = max([max(dof) for dof in self.DOF]) - return nt - - @property - def nm(self): - """No. of kinematic mechanisms, or the no. of dimensions spanned by the - m independent inextensional mechanisms. - - """ - # A = A_matrix(mdl) - pass - - @property - def NOS(self): - nf = self.nf - nq = sum(self.nq) - return nq - nf - - - -class Material(): - def __init__(self, tag, E, nu=None): - self.tag = tag - self.E: float = E - self.elastic_modulus: float = E - self.poisson_ratio = nu - -class XSect(): - def __init__(self, tag, A, I): - self.tag = tag - self.A: float = A - self.I: float = I - -def UnitProperties(): # Legacy; consider removing and adding unit materials to model by default - return (Material('unit', 1.0), XSect('unit', 1.0, 1.0)) - -class State: - """STATE is a data structure with information about the current state of the structure in fields - - """ - def __init__(self, model, method="Linear"): - self.model = model - self.num = len(model.states) - self.data = {"Q": [[qi for qi in elem.q.values()] for elem in model.elems], - "P": {str(dof):0 for dof in [item for sublist in model.DOF for item in sublist]}, - "DOF": 'model.numDOF(model)' - } - self.method = method - - - def eload(self, elem, mag, dirn='y'): - if type(elem) is str: - elem = self.model.delems[elem] - if not(type(mag) is list): - if dirn=='y': - mag = [0.0, mag] - elif dirn=='x': - mag = [mag, 0.0] - elem.w[self.num] = mag - -##>***************************************************************************** - - -# -# ##>***************************************************************************** -# - - -settings = { - "DATAFRAME_LATEX": True, -} - -subscripts = { - 'part': 'subscript', - None: '', - 'initial': '0', - 'continuous':'c', - 'primary':'i', - 'reactions':'d', - 'elem-load':'w', - 'free': 'f', -} - -def del_zeros(mat): - delrows = np.where(~mat.any(axis=1))[0] - delcols = np.where(~mat.any(axis=0))[0] - - newM = np.delete(mat, delcols, axis=1) - newM = np.delete(newM, delrows, axis=0) - return newM - -def _elem_dofs(Elem): - dofs = [] - for node in Elem.nodes: - dofs.extend(node.dofs) - return dofs - -def transfer_vars(item1, item2): - for key, value in item1.__dict__.items(): - item2.__dict__[key] = value - -def Localize(U_vector, P_vector, model=None): - if model is None: - model = U_vector.model - A = A_matrix(model) - Q0 = Q0_vector(model) - Ks = Ks_matrix(model) - - V = A.f @ U_vector.f - Q = Ks@V + Q0 - return V, Q - - - -class Structural_Vector(np.ndarray): - column_data = ["vector"] - row_data = None - subs = [None] - tag = 'Vector' - def __new__(cls, mat): - mat = mat - return np.asarray(mat).view(cls) - - def __init__(self, mat): - self.tag = None - - def _repr_html_(self): - try: - out = self.df.to_html() - return out - except: - return self - - def __add__(self, other): - if isinstance(other, type(self)): - out = np.add(self, other).view(type(self)) - transfer_vars(self, out) - else: - out = super().__add__(other) - return out - - def get(self, key): - idx = np.array([i for i,j in enumerate(self.row_data) if str(j) == key], dtype=int) - # row = self.row_data.index(component) - return self[idx] - - def set_item(self, key, value): - idx = np.array([i for i,j in enumerate(self.row_data) if str(j) == key], dtype=int) - self[idx] = value - - def rows(self, component): - idx = np.where(np.isin(self.row_data,component))[0] - newV = self[idx] - newV.row_data = np.array(self.row_data)[idx] - newV.model = self.model - return newV - - @property - def df(self): - row_data = ['$'+str(tag)+'$' for tag in self.row_data] - header = '$'+self.tag+'_{{' - try: - for sub in self.subs: - header += subscripts[sub] - except: pass - print - header += '}}$' - - return pd.DataFrame(np.around(self,14), index=row_data, columns=[header]) - - @property - def symb(self): - var = [] - for eid in self.row_data: - var.append(sp.symbols(self.tag+eid)) - return sp.Matrix(var) - - @property - def disp(self): - return sp.Matrix(self) - - -class Structural_Matrix(np.ndarray): - column_data = None - row_data = None - c_ridx = None # row indexes for use in .c method - c_cidx = None # column indexes for use in .c method - tag = None - - def __new__(cls, mat): - mat = mat - return np.asarray(mat).view(cls) - - def __init__(self, mat): - self.tag = None - - def __matmul__(self, other): - if isinstance(other, Structural_Matrix): - out = np.matmul(self,other).view(Structural_Matrix) - transfer_vars(self,out) - out.column_data = other.column_data - - elif isinstance(other, Structural_Vector): - out = np.matmul(self,other).view(Structural_Vector) - out.row_data = self.row_data - # out.column_data = ['.'] - else: - out = np.matmul(self,other).view(Structural_Vector) - out.row_data = self.row_data - # out.column_data = ['.'] - return out - - def __add__(self, other): - out = np.add(self,other) - if (isinstance(other, Structural_Matrix) - or isinstance(other, Structural_Vector)): - out.row_data = self.row_data - return out - - - - def __truediv__(self, other): - out = np.ndarray.__truediv__(self, other) - if (isinstance(other, float) - or isinstance(other, Structural_Matrix) - or isinstance(other, Structural_Vector)): - out = np.ndarray.__truediv__(self, other).view(Structural_Matrix) - transfer_vars(self, out) - else: - out = np.ndarray.__truediv__(self, other) - return out - - def _repr_html_(self): - try: - df = self.df - return df.to_html() - except: - try: - return pd.DataFrame(self).to_html() - except: - pass - - @property - def disp(self): - return sp.Matrix(self) - - @property - def df(self): - if settings['DATAFRAME_LATEX']: - row_data = ['$'+str(tag)+'$' for tag in self.row_data] - column_data = ['$'+str(tag)+'$' for tag in self.column_data] - else: - row_data = [str(i) for i in self.row_data] - column_data = [str(i) for i in self.column_data] - return pd.DataFrame(np.around(self,5), index=row_data, columns=column_data) - - @property - def inv(self): - mat = np.linalg.inv(self) - transfer_vars(self, mat) - mat.row_data = self.column_data - mat.column_data = self.row_data - return mat - - @property - def rank(self): - """Return the rank of a matrix""" - return np.linalg.matrix_rank(self) - - @property - def lns(self): - """Return a basis for the left nullspace of a matrix.""" - return scipy.linalg.null_space(self.T) - - @property - def nls(self): - """return a basis for the nullspace of matrix.""" - return scipy.linalg.null_space(self) - - @property - def ker(self): - "Return a basis for the kernel (nullspace) of a matrix." - kernel = scipy.linalg.null_space(self) - ker = Structural_Matrix(kernel) - transfer_vars(self,ker) - ker.row_data = self.column_data - ker.column_data = [str(i+1) for i in range(len(ker[0]))] - return ker - - def lu(self): - - return scipy.linalg.lu(self) - - @property - def c(self): - delcols = self.c_cidx - delrows = self.c_ridx - - newM = np.delete(self, delrows, axis=0).view(type(self)) - if delcols: newM = np.delete(newM, delcols, axis=1).view(type(self)) - - transfer_vars(self, newM) - - if delcols: newM.column_data = list(np.delete(self.column_data, delcols)) - - newM.row_data = list(np.delete(self.row_data, delrows)) - return newM - - def round(self, num): - newM = np.around(self, num).view(Structural_Matrix) - transfer_vars(self,newM) - return newM - - def remove(self, component): - """Remove items by looking up column_data/row_data""" - if type(component) is list: - for item in component: - if item in self.column_data: - delcol = self.column_data.index(item) - newM = np.delete(self, delcol, axis=1).view(type(self)) - transfer_vars(self,newM) - newM.column_data = list(np.delete(self.column_data, delcol)) - else: - delrow = self.row_data.index(item) - newM = np.delete(self, delrow, axis=0).view(type(self)) - transfer_vars(self, newM) - # newM.column_data = self.column_data - # newM.model = self.model - newM.row_data = list(np.delete(self.row_data, delrow)) - - else: - item = component - if item in self.column_data: - delcol = self.column_data.index(item) - newM = np.delete(self, delcol, axis=1).view(type(self)) - newM.row_data = self.row_data - newM.model = self.model - newM.column_data = list(np.delete(self.column_data, delcol)) - # try: newM.rel = self.rel - # except: pass - else: - delrow = self.row_data.index(item) - newM = np.delete(self, delrow, axis=0).view(type(self)) - newM.column_data = self.column_data - newM.model = self.model - newM.row_data = list(np.delete(self.row_data, delrow)) - # try: newM.rel = self.rel - # except: pass - return newM - - def get(self, row_name, col_name): - idxr = np.where(self.row_data == row_name) - idxc = np.where(self.column_data == col_name) - return self[idxr, idxc] - - - def rows(self, component): - rows = [self.row_data.index(item) for item in component] - newM = self[rows,:] - newM.model = self.model - newM.column_data = self.column_data - newM.row_data = list(np.array(self.row_data)[rows]) - return newM - - def del_zeros(self): - """Delete rows and columns of a matrix with all zeros""" - delrows = np.where(~self.any(axis=1))[0] - delcols = np.where(~self.any(axis=0))[0] - - newM = np.delete(self, delcols, axis=1).view(type(self)) - newM = np.delete(newM, delrows, axis=0).view(type(self)) - transfer_vars(self, newM) - - newM.column_data = list(np.delete(self.column_data, delcols)) - newM.row_data = list(np.delete(self.row_data, delrows)) - return newM - - def add_cols(self, component): - if "colinear" in component: - vertical = [elem for elem in self.model.elems if elem.Dx==0.0] - other = set({}) - for elem in self.model.elems: - try: other.add(elem.Dy/elem.Dx) - except ZeroDivisionError: pass - return vertical, other - if type(component) is list: - #ASSUMES component IS A LIST OF COLUMN INDICES - delcols = [self.column_data.index(item) for item in component[1:len(component)]] - i0 = self.column_data.index(component[0]) - newM = np.delete(self, delcols, axis=1).view(type(self)) - for col in delcols: - newM[:,i0] += self[:,col] - newM.row_data = self.row_data - newM.model = self.model - newM.column_data = list(np.delete(self.column_data, delcols)) - # try: newM.rel = self.rel - # except: pass - return newM - - def add_rows(self, component): - if type(component) is list: - delrows = [self.row_data.index(item) for item in component[1:len(component)]] - i0 = self.row_data.index(component[0]) - newA = np.delete(self, delrows, axis=0).view(type(self)) - for row in delrows: - newA[i0,:] += self[row,:] - newA.column_data = self.column_data - newA.model = self.model - newA.row_data = list(np.delete(self.row_data, delrows)) - return newA - -class row_vector (Structural_Vector): - def __new__(cls, Matrix): - V = np.zeros((len(Matrix.row_data))) - return np.asarray(V).view(cls) - - def __init__(self, Matrix): - self.tag = "Y" - self.matrix = Matrix - self.row_data = Matrix.row_data - -class column_vector (Structural_Vector): - - def __new__(cls, Matrix, Vector=None): - V = np.zeros((len(Matrix.column_data))) - return np.asarray(V).view(cls) - - def __init__(self, Matrix, Vector=None): - self.tag = "X" - self.matrix = Matrix - self.row_data = Matrix.column_data - if Vector is not None: - for key in Vector.row_data: - self.set_item(key, Vector.rows([key])) - - -class Static_matrix (Structural_Matrix): - """B_MATRIX static matrix of structural model with 2d/3d truss and 2d frame elements - the function forms the static matrix B for all degrees of freedom and - all basic forces of the structural model specified in data structure MODEL; - the function is currently limited to 2d/3d truss and 2d frame elements - - Parameters - --------------- - - model: ema.Model object - - Partitions - ========================================================================================= - - - B.f : nf x ntq - - - B.c : nf x nq - - - B.fc : nf x nq - - - B.i : ni x nq - - - B.x : nx x nq - - where: - - - ni: number of primary (non-redundant) forces. - - nq: number of total, continuous forces. - - nx: number of redundant forces. - - """ - - ranges = { - 'f': 'free dof rows', - 'i': 'primary force columns', - 'x': 'redundant force columns', - 'c': 'continuous force columns (removes releases)', - 'd': 'reaction force columns', - } - - def __new__(cls, model, matrix=None, rng=None): - fullnq = sum([len(elem.rel) for elem in model.elems]) - B = np.zeros((model.nt,fullnq)) - ci = 0 - for elem in model.elems: - dofs = elem.dofs - bg = elem.bg_matrix() - nq = np.size(bg,1) - for j,dof in enumerate(dofs): - B[int(dof)-1,ci:ci+nq] = bg[j,:] - ci = ci+nq - input_array = B - return np.asarray(input_array).view(cls) - - def __init__(self, model, matrix=None, rng=None): - if rng is None: - self.rng = None - self.model = model - self.row_data = np.array([str(dof) for dof in range(1, model.nt+1)]) - self.column_data = np.array([elem.tag+'_'+key for elem in model.elems for key in elem.rel.keys()]) - if matrix is not None: - fullnq = sum([len(elem.rel) for elem in model.elems]) - self[:,:] = np.zeros((model.nt,fullnq)) - for idxr, rw in enumerate(self.row_data): - if rw in matrix.row_data: - for idxc, cl in enumerate(self.column_data): - if cl in matrix.column_data: - self[idxr,idxc] = matrix.get(rw,cl) - - def __matmul__(self, Vector): - if type(Vector) is nForce_vector: - vect = np.matmul(self, Vector).view(iForce_vector) - vect.row_data = self.row_data - vect.matrix = self - elif type(Vector) is iForce_vector: - vect = np.matmul(self, Vector).view(nForce_vector) - vect.row_data = self.row_data - vect.matrix = self - else: - vect = Structural_Matrix.__matmul__(self, Vector) - return vect - - @property - def f(self): - delrows = [idx for idx, dof in enumerate(self.row_data) if int(dof) > self.model.nf] - - newB = np.delete(self, delrows, axis=0).view(type(self)) - transfer_vars(self, newB) - # newB.model = self.model - newB.row_data = list(np.delete(self.row_data, delrows)) - # newB.column_data = self.column_data - return newB - - @property - def i(self): - """Removes rows of B_matrix corresponding to primary (non-redundant) forces""" - Bf = self.f - idx_i = self.model.idx_i - newB = Bf[:,idx_i] - transfer_vars(Bf, newB) - newB.column_data = Bf.column_data[idx_i] - return newB - - @property - def d(self): - """Removes rows corresponding to free dofs""" - delrows = [idx for idx, dof in enumerate(self.row_data) if int(dof) <= self.model.nf] - newB = np.delete(self, delrows, axis=0).view(type(self)) - transfer_vars(self,newB) - newB.row_data = list(np.delete(self.row_data, delrows)) - return newB - - @property - def c(self): - - Bf = self.f - tags = [elem.tag + "_" + rel for elem in self.model.elems for rel in elem.rel if elem.rel[rel]] - delcols = [Bf.column_data.index(tag) for tag in tags] - newB = np.delete(Bf, delcols, axis=1).view(type(self)) - transfer_vars(Bf, newB) - newB.column_data = list(np.delete(Bf.column_data, delcols)) - return newB - - @property - def o(self): - """Remove columns corresponding to element force releases, then delete zeros""" - Bf = self.f - tags = [elem.tag + "_" + rel for elem in self.model.elems for rel in elem.rel if elem.rel[rel]] - # delcols = [idx for idx, rel in enumerate(self.model.rel) if rel==1] - delcols = [Bf.column_data.index(tag) for tag in tags] - newB = np.delete(Bf, delcols, axis=1).view(type(self)) - transfer_vars(Bf, newB) - newB.column_data = list(np.delete(Bf.column_data, delcols)) - newB = newB.del_zeros() - return newB - - @property - def fc(self): - return self.f.c - - @property - def x(self): - """Removes rows of B_matrix corresponding to primary (non-redundant) forces""" - idx_x = self.model.idx_x - newB = self[:,idx_x] - transfer_vars(self, newB) - newB.column_data = self.column_data[idx_x] - return newB - - @property - def barxi(self): - Bx = self.f.x - Bbarxi = self.bari @ -Bx - Bbarxi.column_data = Bx.column_data - return Bbarxi - - @property - def barx(self): - nQ = len(self.column_data) - nx = len(self.model.redundants) - - Bbarxi = self.barxi - - Bbarx = Structural_Matrix(np.zeros((nQ,nx))) - transfer_vars(self, Bbarx) - Bbarx.column_data = Bbarxi.column_data - Bbarx.row_data = self.column_data - - for idxc, cl in enumerate(Bbarx.column_data): - for idxr, rw in enumerate(Bbarx.row_data): - if rw in Bbarxi.row_data: - Bbarx[idxr,idxc] = Bbarxi.get(rw,cl) - elif cl==rw: - Bbarx[idxr, idxc] = 1. - - return Bbarx - - @property - def bari(self): - return self.i.del_zeros().inv - - @property - def c0(self): - Bf = self.f - tags = [elem.tag + "_" + rel for elem in self.model.elems for rel in elem.rel if elem.rel[rel]] - delcols = [Bf.column_data.index(tag) for tag in tags] - newB = Bf - for col in delcols: - newB[:,col] = [0.]*len(Bf[:,0]) - transfer_vars(Bf, newB) - newB.column_data = list(np.delete(Bf.column_data, delcols)) - return newB - -class nStatic_matrix(Structural_Matrix): - ranges = { - 'f': 'free dof rows', - 'i': 'primary force columns', - 'x': 'redundant force columns', - 'c': 'continuous force columns (removes releases)', - 'd': 'reaction force columns',} - def __new__(cls, arry, model, rcdata): - return np.asarray(arry).view(cls) - - def __init__(self, arry, model, rcdata): - self.model = model - self.row_data = rcdata[0] - self.column_data = rcdata[1] - - def __matmul__(self, Vector): - if type(Vector) is nForce_vector: - vect = np.matmul(self, Vector).view(iForce_vector) - vect.row_data = self.row_data - vect.matrix = self - elif isinstance(Vector, iForce_vector): - vect = np.matmul(self, Vector).view(nForce_vector) - vect.row_data = self.row_data - vect.matrix = self - else: - vect = Structural_Matrix.__matmul__(self, Vector) - return vect - - @property - def f(self): - delrows = [idx for idx, dof in enumerate(self.row_data) if int(dof) > self.model.nf] - - newB = np.delete(self, delrows, axis=0).view(type(self)) - transfer_vars(self, newB) - newB.row_data = np.delete(self.row_data, delrows) - return newB - - @property - def i(self): - """Removes rows of B_matrix corresponding to redundant forces""" - # reducedB = self.f.c.del_zeros() - reducedB = self.f.c - rdts = reducedB.model.redundants - tags = [q.elem.tag + '_'+str(q.nature) for q in rdts] - delcols = [reducedB.column_data.index(tag) for tag in tags] - newB = np.delete(reducedB, delcols, axis=1).view(Static_matrix) - transfer_vars(reducedB, newB) - newB.column_data = list(np.delete(reducedB.column_data, delcols)) - return newB - - @property - def d(self): - """Removes rows corresponding to free dofs""" - delrows = [idx for idx, dof in enumerate(self.row_data) if int(dof) <= self.model.nf] - newB = np.delete(self, delrows, axis=0).view(type(self)) - transfer_vars(self,newB) - newB.row_data = list(np.delete(self.row_data, delrows)) - return newB - - @property - def c(self): - """Removes columns corresponding to element hinges/releases""" - Af = self.f - idx_c = self.model.idx_c - newA = Af[:,idx_c] - transfer_vars(Af, newA) - newA.column_data = Af.column_data[idx_c] - return newA - - @property - def o(self): - """Remove columns corresponding to element force releases, then delete zeros""" - Bf = self.f - tags = [elem.tag + "_" + rel for elem in self.model.elems for rel in elem.rel if elem.rel[rel]] - # delcols = [idx for idx, rel in enumerate(self.model.rel) if rel==1] - delcols = [Bf.column_data.index(tag) for tag in tags] - newB = np.delete(Bf, delcols, axis=1).view(type(self)) - transfer_vars(Bf, newB) - newB.column_data = list(np.delete(Bf.column_data, delcols)) - newB = newB.del_zeros() - return newB - - @property - def fc(self): - return self.f.c - - @property - def x(self): - """Removes rows of B_matrix corresponding to primary (non-redundant) forces - - """ - rdts = self.model.redundants - tags = [q.elem.tag + '_'+str(q.nature) for q in rdts] - cols = [self.column_data.index(tag) for tag in tags] - newB = self[:,cols] - transfer_vars(self, newB) - newB.column_data = np.array([self.column_data[col] for col in cols]) - # newB.row_data = self.row_data - # newB.model = self.model - return newB - - @property - def barxi(self): - Bx = self.f.x - Bbarxi = self.bari @ -Bx - Bbarxi.column_data = Bx.column_data - return Bbarxi - - @property - def barx(self): - nQ = len(self.column_data) - nx = len(self.model.redundants) - - Bbarxi = self.barxi - - Bbarx = Structural_Matrix(np.zeros((nQ,nx))) - transfer_vars(self, Bbarx) - Bbarx.column_data = Bbarxi.column_data - Bbarx.row_data = self.column_data - for idxc, cl in enumerate(Bbarx.column_data): - for idxr, rw in enumerate(Bbarx.row_data): - if rw in Bbarxi.row_data: - Bbarx[idxr,idxc] = Bbarxi.get(rw,cl) - elif cl==rw: - Bbarx[idxr, idxc] = 1. - - return Bbarx - - @property - def bari(self): - return self.i.del_zeros().inv - - @property - def c0(self): - Bf = self.f - tags = [elem.tag + "_" + rel for elem in self.model.elems for rel in elem.rel if elem.rel[rel]] - # delcols = [idx for idx, rel in enumerate(self.model.rel) if rel==1] - delcols = [Bf.column_data.index(tag) for tag in tags] - newB = Bf - for col in delcols: - newB[:,col] = [0.]*len(Bf[:,0]) - transfer_vars(Bf, newB) - newB.column_data = list(np.delete(Bf.column_data, delcols)) - # newB.row_data = self.row_data - # newB.rel = self.rel - return newB - -def B_matrix(model, matrix=None, rng=None): - """Returns a Static_matrix object""" - return Static_matrix(model, matrix, rng) - -def nB_matrix(model): - """Returns a Static_matrix object""" - fullnq = sum([len(elem.rel) for elem in model.elems]) - B = np.zeros((model.nt, fullnq)) - ci = 0 - for elem in model.elems: - eid = elem.dofs - bg = elem.bg_matrix() - nq = np.size(bg,1) - for j,eidi in enumerate(eid): - B[int(eidi)-1,ci:ci+nq] = bg[j,:] - ci = ci+nq - input_array = B - # return np.asarray(input_array).view(nStatic_matrix) - matrix = np.asarray(input_array) - return nStatic_matrix(model, matrix, rng) - -def Bh_matrix(model): - """Returns a Static_matrix object""" - fullnq = sum([len(elem.rel) for elem in model.elems]) - B = np.zeros((model.nt, fullnq)) - ci = 0 - for elem in model.elems: - eid = elem.dofs - bg = elem.bg_matrix(Roption=True) - nq = np.size(bg,1) - for j,eidi in enumerate(eid): - B[int(eidi)-1,ci:ci+nq] = bg[j,:] - ci = ci+nq - - row_data = np.array([str(dof) for dof in range(1, model.nt+1)]) - column_data = np.array([elem.tag+'_'+key - for elem in model.elems - for key in elem.rel.keys()]) - rcdata = (row_data, column_data) - return nStatic_matrix(B, model, rcdata) - -class Kinematic_matrix(Structural_Matrix): - """Class for the kinematic matrix of a structural model with 2d/3d truss and 2d frame elements - the function forms the kinematic matrix A for all degrees of freedom and - all element deformations of the structural model specified in data structure MODEL - the function is currently limited to 2d/3d truss and 2d frame elements - - Returns - --------- - - Kinematic matrix - - """ - - ranges = { - 'f': 'free dof columns', - 'i': 'primary force/deformation rows', - 'x': 'redundant force/deformation rows', - 'd': 'reaction force/deformation rows', - 'c': 'continuous (hingeless) force/deformation rows' - } - - def __new__(cls, model, matrix=None,rng=None): - A = np.zeros((sum(model.nv),model.nt)) - ri = 0 - for elem in model.elems: - eid = elem.dofs - ag = elem.ag() - nv = len(elem.v) - for j, eidi in enumerate(eid): - A[ri:ri+nv, int(eidi)-1] = ag[:,j] - ri = ri+nv - input_array = A - return np.asarray(input_array).view(cls) - - def __init__(self, model, matrix=None,rng=None): - if rng is None: - self.rng = None - self.model = model - self.column_data = np.array([str(dof) for dof in range(1,model.nt+1)]) - self.row_data = np.array([elem.tag+'_'+key for elem in model.elems for key in elem.v.keys()]) - self.basic_deformations = np.array([v for elem in model.elems for v in elem.basic_deformations]) - - self.idx_h = [] - if matrix is not None: - fullnq = sum([len(elem.rel) for elem in model.elems]) - self[:,:] = np.zeros((model.nt,fullnq)) - for idxr, rw in enumerate(self.row_data): - if rw in matrix.row_data: - for idxc, cl in enumerate(self.column_data): - if cl in matrix.column_data: - self[idxr,idxc] = matrix.get(rw,cl) - - def __matmul__(self, Vector): - if isinstance(Vector, Deformation_vector): - # print('a') - vect = np.matmul(self, Vector).view(Displacement_vector) - vect.row_data = self.row_data - vect.matrix = self - elif isinstance(Vector, Displacement_vector): - vect = np.matmul(self, Vector).view(Deformation_vector) - vect.row_data = self.row_data - vect.matrix = self - else: - vect = Structural_Matrix.__matmul__(self, Vector) - - return vect - - - def combine(self, component): - if "colinear" in component: - vertical = [elem for elem in self.model.elems if elem.Dx==0.0] - other = set({}) - for elem in self.model.elems: - try: other.add(elem.Dy/elem.Dx) - except ZeroDivisionError: pass - return vertical, other - - if type(component) is list: - ## TO BE DEPRECATED - #ASSUMES component IS A LIST OF COLUMN INDICES - delcols = [self.column_data.index(item) for item in component[1:len(component)]] - i0 = self.column_data.index(component[0]) - newA = np.delete(self, delcols, axis=1).view(Kinematic_matrix) - for col in delcols: - newA[:,i0] += self[:,col] - newA.row_data = self.row_data - newA.model = self.model - newA.column_data = list(np.delete(self.column_data, delcols)) - # newA.rel = self.rel - return newA - - @property - def f(self): - """Removes columns corresponding to fixed dofs""" - delcols = [idx for idx, dof in enumerate(self.column_data) if int(dof) > self.model.nf] - newA = np.delete(self, delcols, axis=1).view(Kinematic_matrix) - transfer_vars(self, newA) - newA.column_data = list(np.delete(self.column_data, delcols)) - return newA - - @property - def i(self): - """Removes rows corresponding to redundant forces""" - Afc = self.f.c - rdts = self.model.redundants - tags = [q.elem.tag +'_'+ str(q.nature) for q in rdts] - delrows = [Afc.row_data.index(tag) for tag in tags] - newA = np.delete(Afc, delrows, axis=0).view(Kinematic_matrix) - transfer_vars(Afc,newA) - newA.row_data = list(np.delete(Afc.row_data, delrows)) - return newA - - @property - def d(self): - """Removes columns corresponding to free dofs""" - delcols = [idx for idx, dof in enumerate(self.column_data) if int(dof) <= self.model.nf] - newA = np.delete(self, delcols, axis=1).view(Kinematic_matrix) - transfer_vars(self,newA) - newA.column_data = list(np.delete(self.column_data, delcols)) - return newA - - - @property - def c(self): - """Removes rows corresponding to element hinges/releases""" - Af = self.f - idx_c = self.model.idx_c - newA = Af[idx_c,:] - transfer_vars(Af, newA) - newA.row_data = Af.row_data[idx_c] - return newA - - @property - def c0(self): - Af = self.f - n_col = len(Af.T) - tags = [elem.tag + "_" + rel for elem in self.model.elems for rel in elem.rel if elem.rel[rel]] - # delcols = [idx for idx, rel in enumerate(self.model.rel) if rel==1] - delrows = np.where(np.isin(Af.row_data,tags)) - # delrows = [Af.row_data.index(tag) for tag in tags] - newA = Af - for rw in delrows: - newA[rw,:] = [0.]*n_col - transfer_vars(Af, newA) - # newA.row_data = list(np.delete(Af.row_data, delrows)) - return newA - - @property - def o(self): - Af = self.f - n_col = len(Af.T) - tags = [elem.tag + "_" + rel for elem in self.model.elems for rel in elem.rel if elem.rel[rel]] - # delcols = [idx for idx, rel in enumerate(self.model.rel) if rel==1] - delrows = [Af.row_data.index(tag) for tag in tags] - newA = Af - for rw in delrows: - newA[rw,:] = [0.]*n_col - newA = newA.del_zeros() - transfer_vars(Af, newA) - # newA.row_data = list(np.delete(Af.row_data, delrows)) - return newA - - @property - def e(self): - Af = self.f - n_col = len(Af.T) - tags = [elem.tag + "_" + rel for elem in self.model.elems for rel in elem.rel if elem.rel[rel]] - # delcols = [idx for idx, rel in enumerate(self.model.rel) if rel==1] - delrows = [Af.row_data.index(tag) for tag in tags] - newA = Af - for rw in delrows: - newA[rw,:] = [0.]*n_col - newA = newA.del_zeros() - transfer_vars(Af, newA) - # newA.row_data = list(np.delete(Af.row_data, delrows)) - return newA - - -class nKinematic_matrix (Structural_Matrix): - """Class for the kinematic matrix of a structural model with 2d/3d truss and 2d frame elements - the function forms the kinematic matrix A for all degrees of freedom and - all element deformations of the structural model specified in data structure MODEL - the function is currently limited to 2d/3d truss and 2d frame elements - - Returns - --------- - Kinematic matrix - - """ - - ranges = { - 'f': 'free dof columns', - 'i': 'primary force/deformation rows', - 'x': 'redundant force/deformation rows', - 'd': 'reaction force/deformation rows', - 'c': 'continuous (hingeless) force/deformation rows' - } - - def __new__(cls, model, matrix=None,rng=None): - A = np.zeros((sum(model.nv),model.nt)) - ri = 0 - for elem in model.elems: - eid = elem.dofs - ag = elem.ag() - nv = len(elem.v) - for j, eidi in enumerate(eid): - A[ri:ri+nv, int(eidi)-1] = ag[:,j] - ri = ri+nv - input_array = A - return np.asarray(input_array).view(cls) - - def __new__(cls, arry, model, rcdata): - return np.asarray(arry).view(cls) - - def __init__(self, arry, model, rcdata): - self.model = model - self.row_data = rcdata[0] - self.column_data = rcdata[1] - self.basic_deformations = np.array([v for elem in model.elems for v in elem.basic_deformations]) - - self.idx_h = [] - - - def __matmul__(self, Vector): - if isinstance(Vector, Deformation_vector): - # print('a') - vect = np.matmul(self, Vector).view(Displacement_vector) - vect.row_data = self.row_data - vect.matrix = self - elif isinstance(Vector, Displacement_vector): - # print('b') - vect = np.matmul(self, Vector).view(Deformation_vector) - # vect = np.matmul(self, Vector).view(Structural_Vector) - vect.row_data = self.row_data - vect.matrix = self - else: - # print('else') - vect = Structural_Matrix.__matmul__(self, Vector) - return vect - - def combine(self, component): - if "colinear" in component: - vertical = [elem for elem in self.model.elems if elem.Dx==0.0] - other = set({}) - for elem in self.model.elems: - try: other.add(elem.Dy/elem.Dx) - except ZeroDivisionError: pass - return vertical, other - - if type(component) is list: - ## TO BE DEPRECATED - #ASSUMES component IS A LIST OF COLUMN INDICES - delcols = [self.column_data.index(item) for item in component[1:len(component)]] - i0 = self.column_data.index(component[0]) - newA = np.delete(self, delcols, axis=1).view(Kinematic_matrix) - for col in delcols: - newA[:,i0] += self[:,col] - newA.row_data = self.row_data - newA.model = self.model - newA.column_data = list(np.delete(self.column_data, delcols)) - # newA.rel = self.rel - return newA - - @property - def f(self): - """Removes columns corresponding to fixed dofs""" - delcols = [idx for idx, dof in enumerate(self.column_data) if int(dof) > self.model.nf] - newA = np.delete(self, delcols, axis=1).view(Kinematic_matrix) - transfer_vars(self, newA) - newA.column_data = list(np.delete(self.column_data, delcols)) - return newA - - @property - def i(self): - """Removes rows corresponding to redundant forces""" - Afc = self.f.c - rdts = self.model.redundants - tags = [q.elem.tag +'_'+ str(q.nature) for q in rdts] - delrows = [Afc.row_data.index(tag) for tag in tags] - newA = np.delete(Afc, delrows, axis=0).view(Kinematic_matrix) - transfer_vars(Afc,newA) - newA.row_data = list(np.delete(Afc.row_data, delrows)) - return newA - - @property - def d(self): - """Removes columns corresponding to free dofs""" - delcols = [idx for idx, dof in enumerate(self.column_data) if int(dof) <= self.model.nf] - newA = np.delete(self, delcols, axis=1).view(Kinematic_matrix) - transfer_vars(self,newA) - newA.column_data = list(np.delete(self.column_data, delcols)) - return newA - - - @property - def c(self): - """Removes rows corresponding to element hinges/releases""" - Af = self.f - idx_c = self.model.idx_c - newA = Af[idx_c,:] - transfer_vars(Af, newA) - newA.row_data = Af.row_data[idx_c] - return newA - - @property - def c0(self): - Af = self.f - n_col = len(Af.T) - tags = [elem.tag + "_" + rel for elem in self.model.elems for rel in elem.rel if elem.rel[rel]] - # delcols = [idx for idx, rel in enumerate(self.model.rel) if rel==1] - delrows = [Af.row_data.index(tag) for tag in tags] - newA = Af - for rw in delrows: - newA[rw,:] = [0.]*n_col - transfer_vars(Af, newA) - # newA.row_data = list(np.delete(Af.row_data, delrows)) - return newA - - @property - def o(self): - Af = self.f - n_col = len(Af.T) - tags = [elem.tag + "_" + rel for elem in self.model.elems for rel in elem.rel if elem.rel[rel]] - # delcols = [idx for idx, rel in enumerate(self.model.rel) if rel==1] - delrows = [Af.row_data.index(tag) for tag in tags] - newA = Af - for rw in delrows: - newA[rw,:] = [0.]*n_col - newA = newA.del_zeros() - transfer_vars(Af, newA) - # newA.row_data = list(np.delete(Af.row_data, delrows)) - return newA - - @property - def e(self): - Af = self.f - n_col = len(Af.T) - tags = [elem.tag + "_" + rel for elem in self.model.elems for rel in elem.rel if elem.rel[rel]] - # delcols = [idx for idx, rel in enumerate(self.model.rel) if rel==1] - delrows = [Af.row_data.index(tag) for tag in tags] - newA = Af - for rw in delrows: - newA[rw,:] = [0.]*n_col - newA = newA.del_zeros() - transfer_vars(Af, newA) - # newA.row_data = list(np.delete(Af.row_data, delrows)) - return newA - -def A_matrix(Domain, matrix=None): - """Returns a Kinematic_matrix object""" - return Kinematic_matrix(Domain,matrix) - - -class Diag_matrix(Structural_Matrix): - """Block diagonal matrix of element flexibility/stiffness matrices for structural model - - - this class represents the block diagonal matrix of element flexibility or stiffness matrices - for a structural model. - - """ - tag = 'F' - def __new__(cls, arry, rc_data, model): - basic_forces = rc_data - arry = np.asarray(arry).view(cls) - arry.basic_forces = basic_forces - arry.model = model - # arry.column_data = [elem.tag+'_'+key for elem in model.elems for key in elem.rel.keys()] - arry.column_data = arry.row_data = np.array([q.elem.tag+'_'+str(q.number) for q in basic_forces]) - # arry.row_data = [elem.tag+'_'+key for elem in model.elems for key in elem.rel.keys()] - return arry - - def __init__(self, arry, rc_data, model): - """Parameters - ========================================================================================= - model - - """ - self.model = model - - def __matmul__(self, Vector): - if type(Vector) is Deformation_vector: - vect = np.matmul(self, Vector).view(iForce_vector) - vect.row_data = self.row_data - vect.matrix = self - elif isinstance(Vector, iForce_vector): - vect = np.matmul(self, Vector).view(Deformation_vector) - # vect.part = 'continuous' - vect.row_data = self.row_data - vect.matrix = self - - else: - vect = Structural_Matrix.__matmul__(self, Vector) - return vect - - @property - def c(self): - """Removes columns corresponding to element hinges/releases""" - idx_c = self.model.idx_c - newA = self[:,idx_c] - newA = newA[idx_c,:] - transfer_vars(self, newA) - newA.basic_forces = self.basic_forces[idx_c] - newA.column_data = self.column_data[idx_c] - newA.row_data = self.row_data[idx_c] - return newA - -def Fs_matrix(model, Roption=True): - """Returns a Flexibility_matrix object""" - if Roption: - f = np.array([elem.f_matrix(Roption) for elem in model.elems]) - basic_forces = np.array([q for elem in model.elems for q in elem.basic_forces if not q.rel]) - else: - f = np.array([elem.f_matrix(Roption) for elem in model.elems]) - basic_forces = np.array([q for elem in model.elems for q in elem.basic_forces]) - Fs = scipy.linalg.block_diag(*f) - - basic_forces = basic_forces - Fs = Diag_matrix(Fs, basic_forces, model) - return Fs - -def Ks_matrix(model): - """Returns a Flexibility_matrix object""" - k = np.array([elem.k_matrix() for elem in model.elems]) - Ks = scipy.linalg.block_diag(*k) - basic_forces = np.array([q for elem in model.elems for q in elem.basic_forces]) - Ks = Diag_matrix(Ks, basic_forces, model) - return Ks - - -class Stiffness_matrix (Structural_Matrix): - """... - Parameters - ========================================================================================= - model - - ----------------------------------------------------------------------------------------- - """ - tag = 'K' - def __new__(cls, arry, model, Roption=None): - - input_array = np.asarray(arry).view(cls) - input_array.model = model - input_array.column_data = [str(dof) for dof in range(1, model.nt+1)] - input_array.row_data = ['P_{'+str(dof)+'}' for dof in range(1, model.nt+1)] - - return input_array - - def __init__(self, arry, model, Roption=None): - self.subs = [None] - pass - - - def __matmul__(self, Vector): - if type(Vector) is Displacement_vector: - vect = np.matmul(self, Vector).view(nForce_vector) - vect.row_data = self.row_data - vect.matrix = self - elif type(Vector) is nForce_vector: - vect = np.matmul(self, Vector).view(Displacement_vector) - vect.row_data = self.row_data - vect.model = self.model - else: - vect = Structural_Matrix.__matmul__(self, Vector) - return vect - - @property - def f(self): - delrows = [idx for idx, dof in enumerate([str(dof) for dof in range(1, self.model.nt+1)]) if int(dof) > self.model.nf] - - newK = np.delete(self, delrows, axis=0) - newK = np.delete(newK, delrows, axis=1).view(type(self)) - transfer_vars(self, newK) - - newK.row_data = list(np.delete(self.row_data, delrows)) - newK.column_data = list(np.delete(self.column_data, delrows)) - return newK - -def K_matrix(Model): - """Returns a Stiffness_matrix object""" - K = np.zeros((Model.nt, Model.nt)) - for elem in Model.elems: - ke = elem.ke_matrix() - for i, dof in enumerate(elem.dofs): - for j, doff in enumerate(elem.dofs): - K[int(dof)-1, int(doff)-1] += ke[i, j] - - return Stiffness_matrix(K, Model, Roption=None) - - - -def Kt_matrix(Model, State): - """Returns a Stiffness_matrix object""" - K = np.zeros((Model.nt, Model.nt)) - for elem in Model.elems: - kt = elem.kt_matrix(State) - for i, dof in enumerate(elem.dofs): - for j, doff in enumerate(elem.dofs): - K[int(dof)-1, int(doff)-1] += kt[i, j] - - return Stiffness_matrix(K, Model, Roption=None) - - - -class Displacement_vector(column_vector): - tag = 'U' - def __new__(cls, Kinematic_matrix, Vector=None): - U = np.zeros((len(Kinematic_matrix.column_data))) - return np.asarray(U).view(cls) - - def __init__(self, Kinematic_matrix, Vector=None): - self.matrix = Kinematic_matrix - self.row_data = Kinematic_matrix.column_data - self.subs = [None] - if Vector is not None: - for key in Vector.row_data: - if key in self.row_data: - self.set_item(key, Vector.rows([key])) - - @property - def f(self): - """Removes rows corresponding to fixed dofs""" - delrows = [idx for idx, dof in enumerate(self.row_data) if int(dof) > self.model.nf] - newU = np.delete(self, delrows, axis=0).view(Displacement_vector) - newU.row_data = list(np.delete(self.row_data, delrows)) - # newU.matrix = self.matrix - newU.model = self.model - return newU - -class nDisplacement_vector(Structural_Vector): - tag = 'U' - def __new__(cls, arry, model, row_data, Vector=None): - input_array = np.asarray(arry).view(cls) - return input_array - - def __init__(self, arry, model, row_data, Vector=None): - self.subs = [None] - self.model = model - self.row_data = row_data - - if Vector is not None: - for key in Vector.row_data: - if key in self.row_data: - self.set_item(key, Vector.rows([key])) - - @property - def f(self): - """Removes rows corresponding to fixed dofs""" - delrows = [idx for idx, dof in enumerate(self.row_data) if int(dof) > self.model.nf] - newU = np.delete(self, delrows, axis=0).view(nDisplacement_vector) - newU.row_data = list(np.delete(self.row_data, delrows)) - newU.model = self.model - return newU - -def U_vector(model, vector=None): - """Returns a Displacement_vector object""" - U = np.zeros(model.nt) - row_data = [str(dof) for dof in range(1,model.nt+1)] - U = nDisplacement_vector(U, model, row_data) - - if vector is not None: - if len(vector)==len(U): - U[:,0] = vector[:] - else: - for key in vector.row_data: - if key in U.row_data: - U.set_item(key, vector.rows([key])) - return U - - -class iForce_vector(Structural_Vector): - tag = 'Q' - def __new__(cls, arry, model, row_data, Vector=None): - return np.asarray(arry).view(cls) - - def __init__(self, arry, model, row_data, Vector=None): - self.model = model - self.subs = [None] - self.row_data = row_data - if Vector is not None: - for key in Vector.row_data: - if key in self.row_data: - self.set_item(key, Vector.get(key)) - #< - @property - def i(self): - """Removes rows corresponding to redundant forces""" - rdts = self.model.redundants - - tags = [q.elem.tag + str(q.nature) for q in rdts] - delrows = [self.row_data.index(tag) for tag in tags] - newQ = np.delete(self, delrows, axis=0).view(iForce_Vector) - transfer_vars(self, newQ) - newQ.subs.append('primary') - newQ.row_data = list(np.delete(self.row_data, delrows)) - return newQ - - @property - def c(self): - """Remove rows corresponding to element hinges/releases""" - idx_c = self.model.idx_c - newQ = self[idx_c] - transfer_vars(self, newQ) - newQ.row_data = self.row_data[idx_c] - return newQ - - @property - def x(self): - """Remove rows of corresponding to primary forces""" - rdts = self.model.redundants - tags = [q.elem.tag + '_'+str(q.nature) for q in rdts] - rows = [self.row_data.index(tag) for tag in tags] - newV = self[rows] - newV.row_data = [self.row_data[row] for row in rows] - transfer_vars(self, newV) - return newV - -def Q_vector(model, vector=None): - """Returns a iForce_vector object""" - - arry = np.zeros((model.nQ,1)) - row_data = np.array([elem.tag+'_'+key for elem in model.elems for key in elem.rel.keys()]) - Q = iForce_vector(arry, model, row_data) - if vector is not None: - if len(vector)==len(Q): - Q[:,0] = vector[:] - else: - for key in vector.row_data: - if key in Q.row_data: - Q.set_item(key, vector.rows([key])) - return Q - -def Q0_vector(model): - """Returns a vector of initial element forces""" - arry = np.concatenate([elem.q0_vector() for elem in model.elems]) - row_data = [elem.tag+'_'+key for elem in model.elems for key in elem.q.keys()] - return iForce_vector(arry, model, row_data) - -def Qpl_vector(model): - """Returns a vector of element plastic capacities""" - - Qp_pos = [elem.Qp['+'][key] for elem in model.elems for key in elem.Qp['+']] - Qp_neg = [elem.Qp['-'][key] for elem in model.elems for key in elem.Qp['-']] - row_data = [elem.tag+'_'+key for elem in model.elems for key in elem.Qp['-']] - - # del_idx = np.where(~Bf.any(axis=0))[0] - # Qp_pos = np.delete(Qp_pos, del_idx) - # Qp_neg = np.delete(Qp_neg, del_idx) - # row_data = np.delete(row_data, del_idx) - column_data = ['Q_{pl}^+', 'Q_{pl}^-'] - - Qpl = nKinematic_matrix(np.array([Qp_pos, Qp_neg]).T, model, (row_data, column_data)) - return Qpl - -def Qp_vector(model): - """Returns a vector of element plastic capacities""" - B = B_matrix(model) - Bf = B.f - Qp_pos = [elem.Qp['+'][key] for elem in model.elems for key in elem.Qp['+']] - Qp_neg = [elem.Qp['-'][key] for elem in model.elems for key in elem.Qp['-']] - row_data = [elem.tag+'_'+key for elem in model.elems for key in elem.Qp['-']] - - del_idx = np.where(~Bf.any(axis=0))[0] - Qp_pos = np.delete(Qp_pos, del_idx) - Qp_neg = np.delete(Qp_neg, del_idx) - row_data = np.delete(row_data, del_idx) - column_data = ['Q_{pl}^+', 'Q_{pl}^-'] - - Qpl = nStatic_matrix(np.array([Qp_pos, Qp_neg]).T, model, (row_data, column_data)) - return Qpl - -class nForce_vector(Structural_Vector): - tag = 'P' - def __new__(cls, arry, model, row_data, Vector=None): - return np.asarray(arry).view(cls) - - def __init__(self, arry, model, row_data, Vector=None): - self.model = model - self.subs = [None] - self.row_data = row_data - if Vector is not None: - for key in Vector.row_data: - if key in self.row_data: - self.set_item(key, Vector.get(key)) - - - @property - def f(self): - delrows = [idx for idx, dof in enumerate(self.row_data) if int(dof) > self.model.nf] - - newP = np.delete(self, delrows, axis=0).view(type(self)) - transfer_vars(self, newP) - # newB.model = self.model - newP.row_data = list(np.delete(self.row_data, delrows)) - # newB.column_data = self.column_data - return newP - - @property - def d(self): - """Removes rows corresponding to free dofs""" - delrows = [idx for idx, dof in enumerate(self.row_data) if int(dof) <= self.model.nf] - newP = np.delete(self, delrows, axis=0).view(type(self)) - newP.model = self.model - newP.row_data = list(np.delete(self.row_data, delrows)) - newP.subs.append('reactions') - return newP - - #< - -def P_vector(model, vector=None): - P = np.zeros(model.nt) - for node in model.nodes: - p = node.p_vector() - for i, dof in enumerate(node.dofs): - P[int(dof)-1] += p[i] - row_data = [str(dof) for dof in range(1, model.nt+1)] - return nForce_vector(P, model, row_data) - - - -# def P0_vector(model): -# """Returns a _ object""" -# arry = np.concatenate([elem.p0_vector() for elem in model.elems]) -# row_data = [elem.tag+'_'+key for elem in model.elems for key in elem.v.keys()] -# return nForce_vector(arry, model, row_data) - -def P0_vector(model): - P = np.zeros(model.nt) - for elem in model.elems: - dofs = elem.dofs - if hasattr(elem, 'q0_vector'): - p0 = elem.bg_matrix()@elem.q0_vector() - else: - p0 = np.zeros(len(dofs)) - - for i,df in enumerate(dofs): - P[int(df)-1] += p0[i] - - row_data = [str(dof) for dof in range(1, model.nt+1)] - return nForce_vector(P, model, row_data) - -def Pw_vector(model): - P = np.zeros(model.nt) - for elem in model.elems: - dofs = elem.dofs - if len(dofs)==6: - P[int(dofs[0])-1] += elem.w['y']*elem.L/2*elem.sn - P[int(dofs[1])-1] += -elem.w['y']*elem.L/2*elem.cs - P[int(dofs[3])-1] += elem.w['y']*elem.L/2*elem.sn - P[int(dofs[4])-1] += -elem.w['y']*elem.L/2*elem.cs - else: - pw = elem.pw_vector() - for i,df in enumerate(dofs): - P[int(df)-1] += pw[i] - - - row_data = [str(dof) for dof in range(1, model.nt+1)] - return nForce_vector(P, model, row_data) - - -class Deformation_vector(Structural_Vector): - tag = "V" - def __new__(cls, arry, model, row_data, Vector=None): - input_array = np.asarray(arry).view(cls) - return input_array - - def __init__(self, arry, model, row_data, Vector=None): - self.model = model - self.row_data = row_data - - @property - def c(self): - """Removes rows corresponding to element hinges/releases""" - idx_c = self.model.idx_c - newQ = self[idx_c] - transfer_vars(self, newQ) - newQ.row_data = self.row_data[idx_c] - return newQ - - @property - def i(self): - """Removes rows corresponding to redundant forces""" - rdts = self.model.redundants - tags = [q.elem.tag + '_'+str(q.nature) for q in rdts] - delrows = [self.row_data.index(tag) for tag in tags] - newV = np.delete(self, delrows, axis=0).view(type(self)) - transfer_vars(self,newV) - newV.row_data = list(np.delete(self.row_data, delrows)) - newV.subs.append('primary') - return newV - - @property - def x(self): - """Removes rows of corresponding to primary forces - - """ - rdts = self.model.redundants - tags = [q.elem.tag + '_'+str(q.nature) for q in rdts] - rows = [self.row_data.index(tag) for tag in tags] - newV = self[rows] - newV.row_data = [self.row_data[row] for row in rows] - transfer_vars(self, newV) - return newV - -def V_vector(model, vector=None): - """Returns a Deformation_vector object""" - arry = np.zeros((model.nQ,1)) - row_data = np.array([elem.tag+'_'+key for elem in model.elems for key in elem.rel.keys()]) - V = Deformation_vector(arry, model, row_data) - if vector is not None: - if len(vector)==len(V): - V[:,0] = vector[:] - else: - for key in vector.row_data: - if key in V.row_data: - V.set_item(key, vector.rows([key])) - return V - -def V0_vector(model): - """Returns a Deformation_vector object""" - arry = np.concatenate([elem.v0_vector() for elem in model.elems]) - row_data = [elem.tag+'_'+key for elem in model.elems for key in elem.v.keys()] - return Deformation_vector(arry, model, row_data) - - -def Aub_matrix(model, alpha): - """Return the interaction upperbound matrix""" - aub = np.array([elem.aub_matrix(alpha) for elem in model.elems]) - A = scipy.linalg.block_diag(*aub) - return A - diff --git a/ops/opensees/repl/__init__.py b/SRC/opensees/repl/__init__.py similarity index 100% rename from ops/opensees/repl/__init__.py rename to SRC/opensees/repl/__init__.py diff --git a/ops/opensees/repl/cmdshell.py b/SRC/opensees/repl/cmdshell.py similarity index 100% rename from ops/opensees/repl/cmdshell.py rename to SRC/opensees/repl/cmdshell.py diff --git a/ops/opensees/repl/ptkshell.py b/SRC/opensees/repl/ptkshell.py similarity index 100% rename from ops/opensees/repl/ptkshell.py rename to SRC/opensees/repl/ptkshell.py diff --git a/ops/opensees/repl/unix.py b/SRC/opensees/repl/unix.py similarity index 100% rename from ops/opensees/repl/unix.py rename to SRC/opensees/repl/unix.py diff --git a/SRC/opensees/scripts.py b/SRC/opensees/scripts.py deleted file mode 100644 index d40ba1635e..0000000000 --- a/SRC/opensees/scripts.py +++ /dev/null @@ -1,53 +0,0 @@ -# -# Analysis -# -def eigen(script: str, modes=1, verbose=False): - interp = TclInterpreter() - interp.eval(f""" - - {script} - - set options(-verbose) {int(verbose)} - set options(-numModes) {modes} - set options(-file) /dev/stdout - - set PI 3.1415159 - set DOFs {{1 2 3 4 5 6}} - set nodeList [getNodeTags] - - """ + """ - # Initialize variables `omega`, `f` and `T` to - # empty lists. - foreach {omega f T recorders} {{} {} {} {}} {} - - for {set k 1} {$k <= $options(-numModes)} {incr k} { - lappend recorders [recorder Node -node {*}$nodeList -dof {*}$DOFs "eigen $k";] - } - - set eigenvals [eigen $options(-numModes)]; - - set T_scale 1.0 - foreach eig $eigenvals { - lappend omega [expr sqrt($eig)]; - lappend f [expr sqrt($eig)/(2.0*$PI)]; - lappend T [expr $T_scale*(2.0*$PI)/sqrt($eig)]; - } - - # print info to `stdout`. - #if {$options(-verbose)} { - # # puts "Angular frequency (rad/s): $omega\n"; - # # puts "Frequency (Hz): $f\n"; - # # puts "Periods (sec): $T\n"; - #} - - if {$options(-file) != 0} { - source /home/claudio/brace/Scripts/OpenSeesScripts/brace2.tcl - brace2::io::write_modes $options(-file) $options(-numModes) - } - - foreach recorder $recorders { - remove recorder $recorder - } - """) - return interp - diff --git a/ops/opensees/section/__init__.py b/SRC/opensees/section/__init__.py similarity index 100% rename from ops/opensees/section/__init__.py rename to SRC/opensees/section/__init__.py diff --git a/ops/opensees/section/aisc_imperial.py b/SRC/opensees/section/aisc_imperial.py similarity index 100% rename from ops/opensees/section/aisc_imperial.py rename to SRC/opensees/section/aisc_imperial.py diff --git a/ops/opensees/section/analysis.py b/SRC/opensees/section/analysis.py similarity index 100% rename from ops/opensees/section/analysis.py rename to SRC/opensees/section/analysis.py diff --git a/ops/opensees/section/patch.py b/SRC/opensees/section/patch.py similarity index 100% rename from ops/opensees/section/patch.py rename to SRC/opensees/section/patch.py diff --git a/ops/opensees/section/section.py b/SRC/opensees/section/section.py similarity index 100% rename from ops/opensees/section/section.py rename to SRC/opensees/section/section.py diff --git a/ops/opensees/section/torsion.py b/SRC/opensees/section/torsion.py similarity index 100% rename from ops/opensees/section/torsion.py rename to SRC/opensees/section/torsion.py diff --git a/SRC/opensees/state.py b/SRC/opensees/state.py deleted file mode 100644 index 1dec022310..0000000000 --- a/SRC/opensees/state.py +++ /dev/null @@ -1,4 +0,0 @@ -class State: - def __init__(self, interp): - pass - diff --git a/ops/opensees/units/__init__.py b/SRC/opensees/units/__init__.py similarity index 100% rename from ops/opensees/units/__init__.py rename to SRC/opensees/units/__init__.py diff --git a/ops/opensees/units/__main__.py b/SRC/opensees/units/__main__.py similarity index 100% rename from ops/opensees/units/__main__.py rename to SRC/opensees/units/__main__.py diff --git a/ops/opensees/units/cgs.py b/SRC/opensees/units/cgs.py similarity index 100% rename from ops/opensees/units/cgs.py rename to SRC/opensees/units/cgs.py diff --git a/ops/opensees/units/english.py b/SRC/opensees/units/english.py similarity index 100% rename from ops/opensees/units/english.py rename to SRC/opensees/units/english.py diff --git a/ops/opensees/units/fas.py b/SRC/opensees/units/fas.py similarity index 100% rename from ops/opensees/units/fas.py rename to SRC/opensees/units/fas.py diff --git a/ops/opensees/units/fks.py b/SRC/opensees/units/fks.py similarity index 100% rename from ops/opensees/units/fks.py rename to SRC/opensees/units/fks.py diff --git a/ops/opensees/units/fps.py b/SRC/opensees/units/fps.py similarity index 100% rename from ops/opensees/units/fps.py rename to SRC/opensees/units/fps.py diff --git a/ops/opensees/units/ias.py b/SRC/opensees/units/ias.py similarity index 100% rename from ops/opensees/units/ias.py rename to SRC/opensees/units/ias.py diff --git a/ops/opensees/units/iks.py b/SRC/opensees/units/iks.py similarity index 100% rename from ops/opensees/units/iks.py rename to SRC/opensees/units/iks.py diff --git a/ops/opensees/units/ips.py b/SRC/opensees/units/ips.py similarity index 100% rename from ops/opensees/units/ips.py rename to SRC/opensees/units/ips.py diff --git a/ops/opensees/units/mks.py b/SRC/opensees/units/mks.py similarity index 100% rename from ops/opensees/units/mks.py rename to SRC/opensees/units/mks.py diff --git a/ops/opensees/units/si.py b/SRC/opensees/units/si.py similarity index 100% rename from ops/opensees/units/si.py rename to SRC/opensees/units/si.py diff --git a/ops/opensees/units/units.py b/SRC/opensees/units/units.py similarity index 100% rename from ops/opensees/units/units.py rename to SRC/opensees/units/units.py diff --git a/SRC/opensees/validation.py b/SRC/opensees/validation.py deleted file mode 100644 index 16548a10da..0000000000 --- a/SRC/opensees/validation.py +++ /dev/null @@ -1,11 +0,0 @@ - # Command # subcommand # syntax validator -{ - "section": { - "Fiber": ..., - "_fall": ( - (re.compile("fiberSec"), "Fiber") - ) - } -} - - diff --git a/COPYRIGHT b/about/COPYRIGHT similarity index 100% rename from COPYRIGHT rename to about/COPYRIGHT diff --git a/ops/opensees/analysis.py b/ops/opensees/analysis.py deleted file mode 100644 index b6a88bc052..0000000000 --- a/ops/opensees/analysis.py +++ /dev/null @@ -1,134 +0,0 @@ -from . import tcl - -class Analysis: - def __init__(self, model, strategy, patterns, recorders=None): - # - # Runtime - # - if not hasattr(model, "_rt") or model._rt is None: - from . import tcl - self.rt = tcl.TclRuntime() - self.rt.send(model) - else: - self.rt = model - - # - # Strategy - # - # Currently, all analysis constructors take the `strategy` - # argument as a C++ std::vector, so - # parameters must be cast to strings. - self._strategy = { - k: [str(i) for i in v] - for k,v in strategy.items() - } - - # - # Patterns - # - self.patterns = patterns - if patterns is not None: - for pattern in patterns.values(): - self.rt.eval(tcl.dumps(pattern)) - - # - # Recorders - # - if recorders is not None: - self.recorders = { - self.rt.eval(tcl.dumps(recorder)): recorder - # recorder.link(runtime=self.rt): recorder - for recorder in recorders - } - - def system(self, *args): pass - - def integrator(self, *args): pass - - def handler(self, *args): pass - - def analyze(self, *args, **kwds): - self._analysis.analyze(*args) - - - -class StaticAnalysis(Analysis): - def __init__(self, model, strategy=None, patterns=None): - if strategy is None: - strategy = { - "algorithm": ["Newton"], - "integrator": ["LoadControl"] - } - super().__init__(model, strategy=strategy, patterns=patterns) - - # Import C++ bindings and create an instance of the analysis - from . import OpenSeesPyRT as libOpenSeesRT - self._analysis = libOpenSeesRT._StaticAnalysis(self.rt._rt, self._strategy) - - - - - -class DirectIntegrationAnalysis(Analysis): - """ - DirectIntegrationAnalysis(patterns=[ - pattern.MultipleSupport( - components = [ - # node, dof, history - ( 1, 1, pattern.ResponseComponent(accel, displ, veloc)), - ( 2, 1, pattern.ResponseComponent(accel, displ, veloc)) - ] - ), - - pattern.UniformAcceleration([ - # dof history - ( 1, pattern.ResponseComponent(accel, displ, veloc)) - ] - ) - - ]) - """ - - def __init__(self, - model, - patterns: dict = None, - strategy: dict = None, - recorders: list = None, - inherit: str = None, - gravity: dict = None, - time_step: float= None, - steps: int = 1 - ): - if strategy is None: - strategy = { - # "algorithm": ["Newton"], - # "integrator": ["Newmark"] - } - super().__init__(model, strategy=strategy, patterns=patterns, recorders=recorders) - # Import C++ bindings and create an instance of the analysis - from . import OpenSeesPyRT as libOpenSeesRT - self._analysis = libOpenSeesRT._DirectIntegrationAnalysis(self.rt._rt, self._strategy) - - # self.steps, self.time_step = min(( - # (len(pat.series.values), pat.series.time_step) - # for pat in self.patterns.values() - # if pat.series.time_step is not None - # ), - # key = lambda i: i[1] - # ) if time_step is None else (steps, time_step) - - - def analyze(self, steps=None, time_step=None, **kwds): - time_step = time_step or kwds.get("dt", None) or self.time_step - if time_step is None: - raise ValueError("Unable to deduce time step size") - - steps = steps or self.steps - if steps is None: - raise ValueError("Unable to deduce time step size") - while True: - self._analysis.analyze(1, time_step) - else: - return self._analysis.analyze(steps, time_step) - - diff --git a/ops/opensees/openseespy.py b/ops/opensees/openseespy.py deleted file mode 100644 index 2b394e7d7c..0000000000 --- a/ops/opensees/openseespy.py +++ /dev/null @@ -1,790 +0,0 @@ -#===----------------------------------------------------------------------===# -# -# STAIRLab -- STructural Artificial Intelligence Laboratory -# Berkeley, CA -# -#===----------------------------------------------------------------------===# -# -""" -This module implements the OpenSeesPy interface. -Imports can be performed exactly as one would -from openseespy, for example: - ->>> import opensees.openseespy as ops - ->>> from opensees.openseespy import node, model - ->>> from opensees.openseespy import * - -""" -import json -from functools import partial - -from .tcl import Interpreter, _lift - -# something to compare the output of model.analyze to: -successful = 0 - - -class OpenSeesError(Exception): - pass - -def _split_iter(source, sep=None, regex=False): - """ - generator version of str.split() - - :param source: - source string (unicode or bytes) - - :param sep: - separator to split on. - - :param regex: - if True, will treat sep as regular expression. - - :returns: - generator yielding elements of string. - """ - if sep is None: - # mimic default python behavior - source = source.strip() - sep = "\\s+" - if isinstance(source, bytes): - sep = sep.encode("ascii") - regex = True - - if regex: - # version using re.finditer() - if not hasattr(sep, "finditer"): - sep = re.compile(sep) - start = 0 - for m in sep.finditer(source): - idx = m.start() - assert idx >= start - yield source[start:idx] - start = m.end() - yield source[start:] - - else: - # version using str.find(), less overhead than re.finditer() - sepsize = len(sep) - start = 0 - while True: - idx = source.find(sep, start) - if idx == -1: - yield source[start:] - return - yield source[start:idx] - start = idx + sepsize - -class Surface: - def __init__(self, nodes, cells, child, points, split): - import shps.child, shps.plane - self.nodes = nodes - self.cells = cells - self.points = points - self.split = split - self.outline = shps.child.IsoparametricMap(shps.plane.Q9, nodes=points) - - - def walk_edge(self): - import numpy as np - nx, ny = self.split - - nat_exterior = [ - *[ ( x, -1) for x in np.linspace(-1, 1, nx+1)[:-1]], - *[ ( 1, y) for y in np.linspace(-1, 1, ny+1)[:-1]], - *[ ( x, 1) for x in reversed(np.linspace(-1, 1, nx+1)[1:])], - *[ (-1, y) for y in reversed(np.linspace(-1, 1, ny+1)[1:])], - ] - - def find_node(coord): - for tag, xyz in self.nodes.items(): - if np.linalg.norm(np.array(xyz) - np.array(coord)) <= 1e-12: - return tag - - exterior_coords = [self.outline.coord(x) for x in nat_exterior] - - for i in range(1,len(exterior_coords)): - yield (find_node(exterior_coords[i-1]), - find_node(exterior_coords[i])) - - yield (find_node(exterior_coords[-1]), - find_node(exterior_coords[ 0])) - - - def __getitem__(self, item): - pass - - -class OpenSeesPy: - """ - This class is meant to be instantiated as a global singleton - that is private to this Python module. - - It encapsulates an instance of Interpreter which implements an - OpenSees state. - """ - def __init__(self, *args, save=False, echo_file=None, **kwds): - self._interp = Interpreter(*args, **kwds) - self._partial = partial - self._save = save - self._echo = echo_file - - self._mesh = {"line": {}, "quad": {}} - - # Enable OpenSeesPy command behaviors - self._interp.eval("pragma openseespy") - - def _str_call(self, proc_name: str, *args, _final=None, **kwds)->str: - """ - Invoke the Interpreter's eval method, calling - a procedure named `proc_name` with arguments - from args and kwds, after converting Python semantics - to Tcl semantics (via _as_str_arg). - - For example, key-word arguments contained in the `kwds` - dict are converted to a sequence of "-key" and "value" - strings. - """ - - tcl_args = (_as_str_arg(i) for i in args) - tcl_kwds = ( - (f"-{key.replace('_','-')}" if val else "") if isinstance(val, bool) - else f"-{key} " + _as_str_arg(val) - for key, val in kwds.items() - ) - cmd = f"{proc_name} {' '.join(tcl_args)} {' '.join(tcl_kwds)}" - - if _final is not None: - cmd += _as_str_arg(_final) - - # TODO: make sure errors print nicely - try: - ret = self.eval(cmd) - except Exception as e: - raise OpenSeesError() from e - - if ret is None or ret == "": - return None - - parts = ret.split() - # Use json parse to cast return values from string. - # This is faster than the standard ast module. - if len(parts) > 1: - try: return list(map(json.loads, parts)) #json.loads("[" + ",".join(parts) + "]") -# try: return json.loads("[" + ",".join(parts) + "]") - except: return ret - - elif proc_name == "eigen": - # "eigen" should always return a list - return [float(ret)] - - else: - try: return json.loads(ret) - except: return ret - - - def eval(self, cmd: str) -> str: - "Evaluate a Tcl command" - if self._echo is not None: - print(cmd, file=self._echo) - return self._interp.eval(cmd) - - def block3D(self, *args, **kwds): - if isinstance(args[6], list) or isinstance(args[7], dict): - return self._str_call("block3D", *args, **kwds) - - # We have to imitate the OpenSeesPy parser, which - # *requires* hard-coding the number of element args - # expected by each element type. This is terribly - # unstable and limited and should only be used when - # backwards compatibility with the original OpenSeesPy - # is absolutely necessary. - elem_name = args[4] - elem_argc = 7 - elem_args = args[6] - - nl = '\n' - ndm = self._str_call("getNDM") - # loop over remaining args to form node coords - node_args = f"""{{ - {nl.join(" ".join(map(str,args[elem_argc+i*(ndm+1):elem_argc+(i+1)*(ndm+1)])) for i in range(int(len(args[elem_argc:])/(ndm+1))))} - }}""" - - return self._str_call("block3D", *args[:6], elem_args, node_args) - - - def block2D(self, *args, **kwds): - if isinstance(args[5], list): - return self._str_call("block2D", *args, **kwds) - - # We have to imitate the OpenSeesPy parser, which - # *requires* hard-coding the number of element args - # expected by each element type. This is terribly - # unstable and limited and should only be used when - # backwards compatibility with the original OpenSeesPy - # is absolutely necessary. - elem_name = args[4] - elem_argc = { - "quad": 9, - "stdquad": 9, - - "shell": 7, - "shellmitc4": 7, - - "shellnldkgq": 7, - "shelldkgq": 7, - - "bbarquad": 8, - - "enhancedquad": 9, - - "sspquad": 9 - }[elem_name.lower()] -1 - - elem_args = list(args[5:elem_argc]) - - nl = '\n' - ndm = self._str_call("getNDM") - # loop over remaining args to form node coords - node_args = f"""{{ - {nl.join(" ".join(map(str,args[elem_argc+i*(ndm+1):elem_argc+(i+1)*(ndm+1)])) for i in range(int(len(args[elem_argc:])/(ndm+1))))} - }}""" - - return self._str_call("block2D", *args[:5], elem_args, node_args) - - - def timeSeries(self, *args, **kwds): - """ - ['Path', 1, '-values', 0.0, 5.0, 8.0, 7.0, 5.0, 3.0, 2.0, 1.0, 0.0, '-time', 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8] - ['Path', 1, '-values', [0.0, 5.0, 8.0, 7.0, 5.0, 3.0, 2.0, 1.0, 0.0], '-time', [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8]] - """ - - args = list(args) - if "-values" in args: - iv = args.index("-values") - # Count the number of floating-point arguments - for nv, value in enumerate(args[iv+1:]): - if not isinstance(value, float): - nv += 1 - break - else: - # if we didnt break out of the for loop - nv += 2 - - values = args[iv+1:iv+nv] - args = [a for a in args[:iv+1]] + [values] + [a for a in args[iv+nv:]] - - if "-time" in args: - it = args.index("-time") - for nt, value in enumerate(args[it+1:]): - if not isinstance(value, float): - nt += 1 - break - else: - # if we didnt break out of the for loop - nt += 2 - - time = args[it+1:it+nt] - args = [a for a in args[:it+1]] + [time] + [a for a in args[it+nt:]] - - return self._str_call("timeSeries", *args, **kwds) - - def pattern(self, *args, load=None, **kwds): - self._current_pattern = args[1] - - if load is not None: - loads = [ - ("load", k, *v) for k,v in load.items() - ] - return self._str_call("pattern", *args, **kwds, _final=loads) - else: - return self._str_call("pattern", *args, **kwds) - - def load(self, *args, pattern=None, load=None, **kwds): - if pattern is None: - pattern = self._current_pattern - - return self._str_call("nodalLoad", *args, "-pattern", pattern, **kwds) - - def mesh(self, type, tag: int, *args, **kwds): - if type == "line": - return self._mesh_line(tag, 2, args[1:3], *args[3:7], args[7:]) - - - def _mesh_line(self, tag, numnodes, ndtags, id, ndf:int, meshsize, eleType='', eleArgs=()): - import numpy as np - from itertools import count - - ndI, ndJ = ndtags - add_node = partial(self._str_call, "node") - add_element = partial(self._str_call, "element") - - xi = np.array(self._str_call("nodeCoord", ndI)) - xj = np.array(self._str_call("nodeCoord", ndJ)) - - L = np.linalg.norm(xj - xi) - nn = int(L//meshsize) + 1 - - nodes = [None for _ in range(nn)] - nodes[0] = ndI - nodes[nn-1] = ndJ - - node_tags = set(self._str_call("getNodeTags")) - new_node = filter(lambda i: i not in node_tags, count(1)) - elem_tags = set(self._str_call("getEleTags") or []) - new_elem = filter(lambda i: i not in elem_tags, count(1)) - - for i,x in enumerate(np.linspace(xi, xj, nn, endpoint=True)[1:]): - - node_tag = next(new_node) - add_node(node_tag, *x) - - nodes[i+1] = node_tag - - elem_tag = next(new_elem) - - if i < nn and eleType != '' and eleArgs: - add_element(eleType,elem_tag,nodes[i],nodes[i+1],*eleArgs) - - self._mesh["line"][tag] = nodes - - - - def section(self, type: str, sec_tag: int, *args, **kwds): - self._current_section = sec_tag - # TODO: error handling - - if "shape" in kwds: - from opensees.section import from_shape - ndm = int(self.eval("getNDM")) - # kwds["shape"] looks like ("W14X90", matTag, (20,4), units?) - shape = from_shape(type, *kwds.pop("shape"), ndm=ndm) - else: - shape = None - - ret = self._str_call("section", type, sec_tag, *args, **kwds) - - if shape is not None: - for fiber in shape.fibers: - self._str_call("fiber", *fiber.coord, fiber.area, fiber.material, section=sec_tag) - - return ret - - def patch(self, *args, **kwds): - section = self._current_section - return self._str_call("patch", *args, "-section", section, **kwds) - - def layer(self, *args, **kwds): - section = self._current_section - return self._str_call("layer", *args, "-section", section, **kwds) - - def fiber(self, *args, **kwds): - section = self._current_section - return self._str_call("fiber", *args, "-section", section, **kwds) - - - -class Model: - def __init__(self, *args, echo_file=None, **kwds): - self._openseespy = OpenSeesPy(echo_file=echo_file) - self._openseespy._str_call("model", *args, **kwds) - - def export(self, *args, **kwds): - return self._openseespy._interp.export(*args, **kwds) - - def lift(self, type_name: str, tag: int): - return _lift(self._openseespy._interp._tcl.interpaddr(), type_name, tag) - - def asdict(self): - """April 2024""" - return self._openseespy._interp.serialize() - - def setFactor(self, factor): - pass - - def element(self, type, tag, *args, **kwds): - if tag is None: - tag = 1 - ele_tags = self.getEleTags() - if ele_tags is None: - ele_tags = [] - elif isinstance(ele_tags, int): - ele_tags = [ele_tags] - - for existing_tag in ele_tags: - if tag <= existing_tag: - tag = existing_tag + 1 - - return self._openseespy._str_call("element", type, tag, *args, **kwds) - - def getIterationCount(self): - return self._openseespy._str_call("numIter") - - def getResidual(self): - return self._openseespy._str_call("printB", "-ret") - - def getTangent(self, **kwds): - import numpy as np - A = np.array(self._openseespy._str_call("printA", "-ret", **kwds)) - return A.reshape([int(np.sqrt(len(A)))]*2) - - def surface(self, split, element: str=None, args=None, points=None, name=None, kwds=None): - # anchor - # normal - import numpy as np - import shps.plane, shps.block - - add_node = partial(self._openseespy._str_call, "node") - add_element = partial(self._openseespy._str_call, "element") - - if isinstance(element, str): - # element is an element name - cell_type = { - "ShellMITC4": shps.plane.Q4, - }.get(element, shps.plane.Q4) - - elif element is None and name is None: - cell_type = shps.plane.Q4 - element = None - - else: - assert isinstance(name, str) - cell_type = element - element = name - - m_elems = self._openseespy._str_call("getEleTags") - if isinstance(m_elems, int): - m_elems = {m_elems} - - elif m_elems is not None: - m_elems = {tag for tag in m_elems} - - m_nodes = self._openseespy._str_call("getNodeTags") - if isinstance(m_nodes, int): - m_nodes = {m_nodes} - if m_nodes is not None: - m_nodes = { - int(tag): self._openseespy._str_call("nodeCoord", f"{tag}") - for tag in m_nodes - } - - if m_nodes is not None and len(m_nodes) > 0: - join = dict(nodes=m_nodes, cells=m_elems) - else: - join = None - - if kwds is None: - kwds = {} - - nodes, elems = shps.block.block(split, cell_type, points=points, - append=False, join=join, **kwds) - - -# anchor_point, anchor_coord = next(iter(anchor.items())) -# if isinstance(anchor_coord, int): -# anchor_coord = self._openseespy._str_call("nodeCoord", f"{anchor_coord}") - -# anchor_point = np.array([*nodes[anchor_point], 0.0]) - for tag, coord in nodes.items(): - add_node(tag, *coord) #*np.array(coord)-anchor_coord) - - if element is not None: - for tag, elem_nodes in elems.items(): - add_element(element, tag, tuple(map(int,elem_nodes)), *args) - - return Surface(nodes=nodes, cells=elems, child=cell_type, points=points, split=split) - - - def __getattr__(self, name: str): - if name in _OVERWRITTEN: - return getattr(self._openseespy, name) - else: - return self._openseespy._partial(self._openseespy._str_call, name) - - -def _as_str_arg(arg, name: str = None): - """ - Convert arg to a string that represents - Tcl semantics. - """ - import numpy as np - if isinstance(arg, (list,np.ndarray)): - return f"{{{' '.join(_as_str_arg(a) for a in arg)}}}" - - elif isinstance(arg, tuple): - return " ".join(map(str, arg)) - - # parse commands like `section Fiber {...}` - elif isinstance(arg, dict): - return "{\n" + "\n".join([ - f"{cmd} " + " ".join(_as_str_arg(a) for a in val) - for cmd, val in arg.items() - ]) + "}" - - else: - return str(arg) - - - -# The global singleton, for backwards compatibility -_openseespy = OpenSeesPy() - -# A list of symbol names that are importable -# from this module. All of these are dynamically -# resolved by the function __getattr__ below. -__all__ = [ -# - "tcl", - "OpenSeesError", - -# OpenSeesPy attributes - - "uniaxialMaterial", - "testUniaxialMaterial", - "setStrain", - "getStrain", - "getStress", - "getTangent", - "getDampTangent", - "wipe", - "model", - "node", - "fix", - "element", - "timeSeries", - "pattern", - "load", - "system", - "numberer", - "constraints", - "integrator", - "algorithm", - "analysis", - "analyze", - "test", - "section", - "fiber", - "patch", - "layer", - "geomTransf", - "beamIntegration", - "loadConst", - "eleLoad", - "reactions", - "nodeReaction", - "eigen", - "modalProperties", - "responseSpectrumAnalysis", - "nDMaterial", - "block2D", - "block3D", - "rayleigh", - "wipeAnalysis", - "setTime", - "remove", - "mass", - "equalDOF", - "nodeEigenvector", - "getTime", - "setCreep", - "eleResponse", - "sp", - "fixX", - "fixY", - "fixZ", - "reset", - "initialize", - "getLoadFactor", - "build", - "printModel", - "printA", - "printB", - "printGID", - "testNorm", - "testNorms", - "testIter", - "recorder", - "database", - "save", - "restore", - "eleForce", - "eleDynamicalForce", - "nodeUnbalance", - "nodeDisp", - "setNodeDisp", - "nodeVel", - "setNodeVel", - "nodeAccel", - "setNodeAccel", - "nodeResponse", - "nodeCoord", - "setNodeCoord", - "getPatterns", - "getFixedNodes", - "getFixedDOFs", - "getConstrainedNodes", - "getConstrainedDOFs", - "getRetainedNodes", - "getRetainedDOFs", - "updateElementDomain", - "getNDM", - "getNDF", - "eleNodes", - "eleType", - "nodeDOFs", - "nodeMass", - "nodePressure", - "setNodePressure", - "nodeBounds", - "start", - "stop", - "modalDamping", - "modalDampingQ", - "setElementRayleighDampingFactors", - "region", - "setPrecision", - "searchPeerNGA", - "domainChange", - "record", - "metaData", - "defaultUnits", - "stripXML", - "convertBinaryToText", - "convertTextToBinary", - "getEleTags", - "getCrdTransfTags", - "getNodeTags", - "getParamTags", - "getParamValue", - "sectionForce", - "sectionDeformation", - "sectionStiffness", - "sectionFlexibility", - "sectionLocation", - "sectionWeight", - "sectionTag", - "sectionDisplacement", - "cbdiDisplacement", - "basicDeformation", - "basicForce", - "basicStiffness", - "InitialStateAnalysis", - "totalCPU", - "solveCPU", - "accelCPU", - "numFact", - "numIter", - "systemSize", - "version", - "setMaxOpenFiles", - "limitCurve", - "imposedMotion", - "imposedSupportMotion", - "groundMotion", - "equalDOF_Mixed", - "rigidLink", - "rigidDiaphragm", - "ShallowFoundationGen", - "setElementRayleighFactors", - "mesh", - "remesh", - "parameter", - "addToParameter", - "updateParameter", - "setParameter", - "getPID", - "getNP", - "barrier", - "send", - "recv", - "Bcast", - "frictionModel", - "computeGradients", - "sensitivityAlgorithm", - "sensNodeDisp", - "sensNodeVel", - "sensNodeAccel", - "sensLambda", - "sensSectionForce", - "sensNodePressure", - "getNumElements", - "getEleClassTags", - "getEleLoadClassTags", - "getEleLoadTags", - "getEleLoadData", - "getNodeLoadTags", - "getNodeLoadData", - "randomVariable", - "getRVTags", - "getRVParamTag", - "getRVValue", - "getMean", - "getStdv", - "getPDF", - "getCDF", - "getInverseCDF", - "correlate", - "performanceFunction", - "gradPerformanceFunction", - "transformUtoX", - "wipeReliability", - "updateMaterialStage", - "sdfResponse", - "probabilityTransformation", - "startPoint", - "randomNumberGenerator", - "reliabilityConvergenceCheck", - "searchDirection", - "meritFunctionCheck", - "stepSizeRule", - "rootFinding", - "functionEvaluator", - "gradientEvaluator", - "getNumThreads", - "setNumThreads", - "logFile", - "setStartNodeTag", - "hystereticBackbone", - "stiffnessDegradation", - "strengthDegradation", - "strengthControl", - "unloadingRule", - "partition", - "pressureConstraint", - "domainCommitTag", -# "runFOSMAnalysis", - "findDesignPoint", - "runFORMAnalysis", - "getLSFTags", - "runImportanceSamplingAnalysis", - "IGA", - "NDTest", -] - -_PROTOTYPES = { -} - -# Commands that are pre-processed in Python -# before forwarding to the Tcl interpreter -_OVERWRITTEN = { - "timeSeries", - "pattern", "load", - "eval", - "section", "patch", "layer", "fiber", - "block2D", - "block3D", - "mesh" -} - - - -def __getattr__(name: str): - # For reference: - # https://peps.python.org/pep-0562/#id4 - if name in _OVERWRITTEN: - return getattr(_openseespy, name) - else: - return _openseespy._partial(_openseespy._str_call, name) - diff --git a/ops/opensees/recorder.py b/ops/opensees/recorder.py deleted file mode 100644 index abb0253859..0000000000 --- a/ops/opensees/recorder.py +++ /dev/null @@ -1,177 +0,0 @@ -from opensees.library.ast import * -from opensees.library.obj import LibCmd - -recorder = LibCmd("recorder") - -class TimeSeries: pass -class Node: pass -class DOF: pass - -node_responses = Str("recorder", about="a string indicating response required. Response types are given in table below.", enum={ - "disp": "displacement*", - "vel": "velocity*", - "accel": "acceleration*", - "incrDisp": "incremental displacement", - "eigen {i}": "eigenvector for mode i", - "reaction": "nodal reaction", - "rayleighForces": "damping forces", - } - ) - -common = [ - Str("destination"), -# Alt("destination", [ -# Str("dest_txt", flag="-file"), -# Str("dest_xml", flag="-xml"), -# Str("dest_bin", flag="-binary"), -# Grp("dest_tcp", flag="-tcp", args=[ -# Str("inetAddr", about='ip address, "xx.xx.xx.xx", of remote machine to which data is sent'), -# Str("port", about='port on remote machine awaiting tcp'), -# ] -# ) -# ], -# reqd=True, -# about="name of file to which output is sent."\ -# "file output is either in xml format (`-xml` option), textual (`-file` option) or binary (`-binary` option)" -# ), - - Int("precision", reqd=False, flag="-precision", - about="number of significant digits (default is 6)"\ - "(optional, default: records at every time step)" - ), - - Flg("-time", reqd=False, - about="using this option places domain time in first entry of each data line, default is to have time omitted" - ), - -# <-closeOnWrite> - Flg("-closeOnWrite", reqd=False, about=""" - using this option will instruct the recorder to invoke a - close on the data handler after every timestep. If this is a file - it will close the file on every step and then re-open it for the - next step. Note, this greatly slows the execution time, but is - useful if you need to monitor the data during the analysis."""), - - Num("time_step", flag="-dT", reqd=False, - about="time interval for recording. will record when next step is `deltaT` greater than last recorder step." - ), -# <-region $regionTag> -# $regionTag a region tag; to specify all nodes in the previously defined region. (optional) - -] - - -class Recorder: - def init(self, format=None): - format = format or self.kwds.get("format", None) - if format is None: - format = self.destination.split(".")[-1] - - if format not in ["txt", "bin", "xml", "binary", "tcp"]: - raise ValueError("Unable to deduce format") - - format = {"txt": "file", "bin": "binary"}.get(format, format) - - self._args[0].flag = "-" + format - - self._runtime = None - - def link(self, *args, runtime): - self._runtime = runtime - return runtime.send(self) - - def parse(self): - if self.dest_txt: - with open(self.dest_txt, "r") as f: - return self.parse_txt(f) - - def update(self): - self._data = self.parse() - - def __getitem__(self, indices): - return self._data[self.index(indices)] - - - -#recorder Node -@recorder -class Node(Recorder): - """ - Record the response of a number of nodes at every converged step. - - - In case you want to remove a recorder, you need to know the tag for that - recorder. Here is an example on how to get the tag of a recorder: - - set tagRc [recorder Node -file nodesD.out -time -node 1 2 3 4 -dof 1 2 disp] - """ - - def parse_xml(self, stream): - pass - - def parse_txt(self, stream): - import numpy as np - return np.loadtxt(stream).reshape((-1, len(self.nodes), len(self.dofs))) - - def parse_bin(self, stream): - pass - - def _index(self, indices): - pass - - @property - def node(self): - nn = len(self.nodes) - ndf = len(self.dofs) - return np.moveaxis(self._data, 1, 0) - - _pysn = [0, -1] - _args = [ - *common, - - Ref("series", flag="-timeSeries", type=TimeSeries, reqd=False, - about="the tag of a previously constructed `TimeSeries`, results from node at each time step are added to load factor from series" - ), - -# <-node $node1 $node2...> - Grp("nodes", type=Ref(type="node"), reqd=False, min=1, flag="-node", about="tags of nodes whose response is being recorded (optional, default: omitted)"), - -# <-nodeRange $startNode $endNode> -# $startNode $endNode.. tag for start and end nodes whose response is being recorded (optional, default: omitted) - -# -dof ($dof1 $dof2 ...) - Grp("dofs", type=Int, reqd=False, num="*", flag="-dof", about="the specified dof at the nodes whose response is requested."), - - node_responses - ] - - - -@recorder -class Element: - _args = [ - *common - ] - -# recorder Element -# """ -# <-file $fileName> <-xml $fileName> <-binary $fileName> -# $fileName name of file to which output is sent. -# file output is either in xml format (-xml option), textual (-file option) or binary (-binary option) -# <-precision $nSD> -# $nSD number of significant digits (optional, default is 6) -# <-time> -# -time (optional using this option places domain time in first entry of each data line, default is to have time ommitted) -# <-closeOnWrite> -# -closeOnWrite optional. using this option will instruct the recorder to invoke a close on the data handler after every timestep. If this is a file it will close the file on every step and then re-open it for the next step. Note, this greatly slows the execution time, but is useful if you need to monitor the data during the analysis. -# <-dT $deltaT> -# $deltaT time interval for recording. will record when next step is $deltaT greater than last recorder step. (optional, default: records at every time step) -# <-ele ($ele1 $ele2 ...)> -# $ele1 $ele2 .. tags of elements whose response is being recorded -- selected elements in domain (optional, default: omitted) -# <-eleRange $startEle $endEle> -# $startEle $endEle .. tag for start and end elements whose response is being recorded -- range of selected elements in domain (optional, default: omitted) -# <-region $regTag> -# $regTag previously-defined tag of region of elements whose response is being recorded -- region of elements in domain (optional) -# $arg1 $arg2 ... -# $arg1 $arg2 ... arguments which are passed to the `setResponse()` element method -# - diff --git a/ops/opensees/series.py b/ops/opensees/series.py deleted file mode 100644 index b86f5aec7a..0000000000 --- a/ops/opensees/series.py +++ /dev/null @@ -1,34 +0,0 @@ - - -import math - -def triangle(t, period=2*math.pi, amplitude=2.0): - a, p = amplitude, period - return 4*a/p * abs((((t-p/4)%p)+p)%p - p/2) - a - - -def sawtooth(t, period=2*math.pi, amplitude = 2): - return amplitude*(t/period - math.floor(0.5 + t/period)) - -def square(t, period=2*math.pi, amplitude=2): - pass - -def stairs(): - pass - -def ramp(): - pass - -def linspace(): - pass - -def points(s_ref, steps): - if isinstance(steps, int): - steps = [steps]*(len(s_ref) - 1) - import numpy as np - diff = np.diff(s_ref,axis=0) - s = np.array([n/stepi*diff[i]+s_ref[i] for i,stepi in enumerate(steps) for n in range(1,stepi+1)]) - return s - -#def cycles(number, fidelity=2, shape="triangle"): -# return getattr(__module__ diff --git a/ops/opensees/tcl.py b/ops/opensees/tcl.py deleted file mode 100644 index 3463fab799..0000000000 --- a/ops/opensees/tcl.py +++ /dev/null @@ -1,390 +0,0 @@ -#===----------------------------------------------------------------------===# -# -# STAIRLab -- STructural Artificial Intelligence Laboratory -# Berkeley, CA -# -#===----------------------------------------------------------------------===# -# -import os -import sys -import json -import pathlib -import platform -from contextlib import contextmanager - -try: - # On certain servers (Heroku), tkinter is not - # provided. In this case we can fall pack - # to the tcinter repackaging - import tkinter -except: - import tcinter as tkinter - -from opensees.library.obj import Component - -class TclError(Exception): - pass - -def exec(script: str, silent=False, analysis=True)->dict: - """ - Execute a Tcl script and return interpreter - """ - if silent: - puts = "proc puts {args} {}" - else: - puts = "" - - - interp = Interpreter(safe=True) - - # Turn off analysis - if analysis is False: - interp.eval(f""" - pragma analysis off - proc eigen {{args}} {{set a 1}} - # ANALYZE SHOULD RETURN A NUMBER IN CASE ITS CHECKED IN SCRIPT - proc analyze {{args}} {{set a 1}} - {puts} - """) - - - # Evaluate the script - _script = f""" - # Wrap passed script in procedure so if there is a file-level - # return, it doesnt skip our print command. - proc _sees_script {{}} {{ - {script} - }} - _sees_script - """ - - interp.eval(script) - - return interp - -def eval(script: str): - import textwrap - interp = _create_interp() - return interp.eval(textwrap.dedent(f""" - {script} - """)) - -def dumps(obj, skip_int_refs=False)->str: - -# TODO: Move this function, maybe to emit - - if not isinstance(obj, (Component,list,tuple)): - # Build out a model - from opensees.emit import OpenSeesWriter - return OpenSeesWriter(obj).dump() - - else: - from opensees.emit.opensees import TclScriptBuilder - writer = TclScriptBuilder(skip_int_refs=skip_int_refs) - try: - writer.send(obj) - if not writer.python_objects: - return writer.getScript(indexed=True) - else: - return writer - except Exception as e: - raise e - # raise ValueError("Cannot dump model with binary objects") - -def _lift(interpaddr: int, type, tag: int): - type = type.lower() - # libOpenSeesRT must be imported by Python - # AFTER if has been loaded by Tcl (this was done - # when a TclRuntime() is created) so that Tcl stubs - # are initialized. Otherwise there will be a segfault - # when a python c-binding attempts to call a Tcl - # C function. Users should never import OpenSeesPyRT - # themselves - from opensees import OpenSeesPyRT as libOpenSeesRT - - _builder = libOpenSeesRT.get_builder(interpaddr) - - if type == "uniaxialmaterial": - return _builder.getUniaxialMaterial(int(tag)) - - elif type == "section": - return _builder.getSection(int(tag)) - - elif type == "backbone": - return _builder.getHystereticBackbone(int(tag)) - - else: - raise TypeError(f"Unimplemented type {type}") - - -class Interpreter: - def __init__(self, model=None, verbose=False, safe=False, preload=True, enable_tk=False): - self._tcl = _create_interp(verbose=verbose, - preload=preload, - enable_tk=enable_tk) - - # TODO: - if not safe: - self._tcl.createcommand("=", self._pyexpr) - - if enable_tk: - self._tcl.createcommand("tkloop", self._tcl.mainloop) - - self._tcl.createcommand("export", self.export) - - # self._tcl.createcommand("import", self._pyimport) - - try: - import sees.tcl - sees.tcl.add_commands(self) - except: - pass - - if model is not None: - self.send(model) - - def __del__(self): - try: - # _tcl may already be deleted - self.eval("wipe") - except: - pass - - def eval(self, string): - try: - return self._tcl.tk.eval(string) - - except tkinter._tkinter.TclError as e: - raise TclError(self._tcl.getvar("errorInfo")) -# print(self._tcl.getvar("errorInfo"), file=sys.stderr) -# raise e - - def serialize(self)->dict: - import tempfile - tmp = tempfile.NamedTemporaryFile(delete=False) - tmp.close() - - self.eval(f"print -json -file {tmp.name}") - - # Read in the generated JSON - with open(tmp.name, "r") as f: - model = json.load(f) - - # Delete the temporary file - os.remove(tmp.name) - return model - - def export(self, *args): - import io - import opensees.emit.mesh - - model = self.serialize() - - if args[0] == "stdout": - file = io.StringIO() - else: - file = args[0] - - if len(args) > 1: - fmt = args[1] - else: - fmt = "vtk" - - # create a meshio object - mesh = opensees.emit.mesh.dump(model, args[0], fmt) - - try: - mesh.write(file, file_format=fmt) - - except Exception as e: - print(e, file=sys.stderr) - - return "" - - def _pyimport(self, *args): - try: - lib = __import__(args[0]) - print(lib) - except: - self.eval("opensees::import " + " ".join(args)) - return - - def _pyexpr(self, *args): - try: - import numpy as math - except: - import math - - env = math.__dict__ - env["locals"] = None - env["globals"] = None - env["__name__"] = None - env["__file__"] = None - env["__builtins__"] = {} - for k in self.eval("info globals").split(): - try: - update = {k: float(self._tcl.getvar(k))} - except: - continue - env.update(update) - try: - return __builtins__["eval"]((" ".join(args[:])).replace("$",""), env) - - except Exception as e: - print(e, file=sys.stderr) - - - -class ModelRuntime: - def __init__(self, ndm=None, ndf=None, **kwds): - self._interp = Interpreter(**kwds) - self._tcl = self._interp._tcl - self._c_domain = None - self._c_rt = None - if ndm is not None: - self.model(ndm, ndf) - - def model(self, ndm, ndf, **kwds): - """ - model(model: opensees.model) - model(ndm:int, ndf:int) - """ - self._interp.eval(f"model basic -ndm {ndm} -ndf {ndf}") - - def eval(self, cmd: str): - return self._interp.eval(cmd) - - def lift(self, type: str, tag: int=None): - """ - Experimental - """ - if type == "uniaxialmaterial": - self.model(2,3) - return _lift(self._tcl.interpaddr(), type, tag) - - def send(self, obj, ndm=2, ndf=3, **kwds): - self.model(ndm=ndm, ndf=ndf) - - m = dumps(obj) - - if isinstance(m, str): - try: - self._interp.eval(m) - except Exception as e: - raise e - print(e, file=sys.stderr) - else: - self.eval(m.getIndex()) - from . import OpenSeesPyRT as libOpenSeesRT - _builder = libOpenSeesRT.get_builder(self._tcl.interpaddr()) - for ident,obj in m.python_objects.items(): - tag = self._interp.eval(f"set {ident.tclstr()}") - _builder.addPythonObject(int(tag), obj) - - script = m.getScript() - self._interp.eval(script) - - return self - - @property - def _rt(self): - # if self._c_rt is None: - # from . import OpenSeesPyRT as libOpenSeesRT - # self._c_rt = libOpenSeesRT.getRuntime(self._tcl.tk.interpaddr()) - return self._c_rt - - @property - def _domain(self): - if self._c_domain is None: - from . import OpenSeesPyRT as libOpenSeesRT - self._c_domain = libOpenSeesRT.get_domain(self._rt) - return self._c_domain - - def getNodeResponse(self, node, typ): - import numpy as np - return np.array(self._domain.getNodeResponse(node, typ)) - - def getTime(self): - return self._domain.getTime() - - time = getTime - - def fix(self, nodes, *dofs): - if not isinstance(nodes,list): - nodes = [nodes] - for node in nodes: - self._interp.eval(f"fix {node} {' '.join(map(str,dofs))}") - - def add_tagged(self, objs): - for k,v in objs.items(): - if isinstance(k, int): - self._interp.eval(v.cmd) - - -@contextmanager -def _build_extension_env(): - cookies = [] - - # Windows and Python >= 3.8 - if hasattr(os, "add_dll_directory"): - for path in os.environ.get("PATH", "").split(os.pathsep): - - if path and pathlib.Path(path).is_absolute() and pathlib.Path(path).is_dir(): - cookies.append(os.add_dll_directory(path)) - try: - yield - - finally: - for cookie in cookies: - cookie.close() - -def _find_openseesrt(): - if "OPENSEESRT_LIB" in os.environ and len(os.environ["OPENSEESRT_LIB"]): - libOpenSeesRT_path = pathlib.Path(os.environ["OPENSEESRT_LIB"]) - return libOpenSeesRT_path.parents[0], libOpenSeesRT_path - - op_sys = platform.system() - ext = { - "Darwin": ".dylib", - "Linux": ".so", - "Windows": ".dll" - }[op_sys] - - install_dir = pathlib.Path(__file__).parents[0] - - if op_sys == "Windows": - libOpenSeesRT_path = install_dir/"bin"/f"OpenSeesRT{ext}" - else: - libOpenSeesRT_path = install_dir/f"libOpenSeesRT{ext}" - - return libOpenSeesRT_path.parents[0], libOpenSeesRT_path - - -def _create_interp(verbose=False, tcl_lib=None, preload=True, enable_tk=False): - - install_dir, libOpenSeesRT_path = _find_openseesrt() - - if verbose: - print(f"OpenSeesRT: {libOpenSeesRT_path}", file=sys.stderr) - - if enable_tk: - interp = tkinter.Tk() - - else: - interp = tkinter.Tcl() - - - if preload: - with _build_extension_env(): - interp.eval(f"load {{{libOpenSeesRT_path}}}") - - def check(): - interp.after(50, check) - - interp.after(50, check) - - return interp - - -TclRuntime = ModelRuntime; # Interpreter - diff --git a/.clang-format b/tools/.clang-format similarity index 100% rename from .clang-format rename to tools/.clang-format diff --git a/tools/SCRIPTS/strain.tcl b/tools/SCRIPTS/strain.tcl new file mode 100755 index 0000000000..479a1ab47d --- /dev/null +++ b/tools/SCRIPTS/strain.tcl @@ -0,0 +1,389 @@ +#!/usr/bin/env -S python -m opensees --enable-tk + +# OpenSees -- Open System for Earthquake Engineering Simulation +# Pacific Earthquake Engineering Research Center +# +# http://opensees.berkeley.edu/ +# +# Uniaxial Material Viewer +# ------------------------ +# +# To run: +# +# python -m opensees --enable-tk strain.tcl +# +# Written: fmk/MHS/cmp +# Date: March 2001 + +# ############################################################## +# Define some global variables & set their initial values +# ############################################################## + +# current material id tag +set matID 0 + +# initial size of canvas used to draw stress-strain relations +set height 400 +set width 400 +set interpreterWidth 60 +set interpreterHeight 8 + +# absolute max values for strain & stress drawn +set maxStress 150.0 +set maxStrain 0.05 + +# last point plotted on canvas +set xLast [expr $height/2] +set yLast [expr $width/2] + +set minE -$maxStrain +set maxE $maxStrain + +# some frames that will be toggled on & off +set toggleFrame .interpreter + +# ############################################################## +# Set main window parameters +# ############################################################## + +wm title . "Uniaxial Material" + +# ############################################################## +# Invoke OpenSees command to make avaliable commands for testing +# the materials: +# 1)uniaxialMaterial +# 2)strainUniaxialTest +# 3)stressUniaxialTest +# 4)tangUniaxialTest +# ############################################################## + +model basic -ndm 1 + +# ############################################################## +# Define menu frame at top of window +# ############################################################## + +frame .mbar -borderwidth 1 -relief raised +pack .mbar -fill x + +menubutton .mbar.materials -text "Materials" -relief raised -menu .mbar.materials.menu +pack .mbar.materials -side left +button .mbar.values -text Values -command "SetValues" +pack .mbar.values -side left +button .mbar.settings -text Settings -command "Settings" +pack .mbar.settings -side left +button .mbar.reset -text Reset -command "Reset" +pack .mbar.reset -side left +button .mbar.quit -text Quit -command exit +pack .mbar.quit -side left + +frame .style -borderwidth 1 -relief sunken +pack .style -fill x -pady 1 + + +# ############################################################## +# proc for buttons +# ############################################################## + + +set buttonlist { .mbar.materials .mbar.values .mbar.settings .mbar.reset } + +# Procedure Settings - used to disable/enable all the buttons in the buttonlist + +proc disable_buttons {} { + global buttonlist + foreach x $buttonlist { + $x configure -state disabled + } +} + +proc enable_buttons {} { + global buttonlist + foreach x $buttonlist { + $x configure -state normal + } +} + + + +# ############################################################## +# Material Menu +# ############################################################## + +set m [menu .mbar.materials.menu] + +# source in the material data structures and procedures +source uniaxialMaterials.tcl +#source x00_uniaxialMaterials.tcl + +# ############################################################## +# Define sample input data +# ############################################################## + + +proc uniaxialMaterialSample {matName args} { + + disable_buttons + + global matID + incr matID + + set w [toplevel .sampleInput] + wm title ${w} ${matName} + + set count 0 + foreach arg $args { + + set posEq [string first "=" ${arg}] + + global matArg${count} + + if {${posEq} != -1} { + set argName [string range ${arg} 0 [expr ${posEq} - 1]] + set argVal [string range ${arg} [expr ${posEq} + 1] end] + } else { + set argName "" + set argVal ${arg} + } + + label ${w}.l${count} -text ${argName} + entry ${w}.e${count} -textvariable matArg${count} -relief sunken + set matArg${count} ${argVal} + + if {${argName} == "matTag"} { + set matArg${count} ${matID} + ${w}.e${count} configure -state disabled + } elseif {${posEq} == -1} { + ${w}.e${count} configure -state disabled + } + + grid ${w}.l${count} -row ${count} -column 1 -sticky e -padx 5 + grid ${w}.e${count} -row ${count} -column 2 + + incr count + } + + button ${w}.ok -text "OK" -command "doneUniaxialMaterial ${matName} ${count}; destroy ${w}; enable_buttons" -padx 10 + grid ${w}.ok -row ${count} -column 1 -columnspan 2 + +} + +proc doneUniaxialMaterial {matName count} { + + set matCommand "uniaxialMaterial ${matName}" + + for {set i 0} {${i} < ${count}} {incr i} { + set tmp "matArg${i}" + global ${tmp} + eval "set matCommand \"${matCommand} $${tmp}\"" + } + + global matID + set fileOpen [open results.out w] + close $fileOpen + puts "$matCommand" + eval "${matCommand}" + invoke UniaxialMaterial $matID { + # eval uniaxialTest $matID + SetValues + Reset + } + +} + +# ############################################################## +# Define the canvas and slider +# ############################################################## + +frame .figure + +# Create the scale (slider) +set theScale [ scale .figure.scale -from [expr -$maxStrain] -to $maxStrain \ + -length $width -variable strain -orient horizontal \ + -tickinterval [expr $maxStrain / 2] -resolution [expr $maxStrain / 100] \ + -showvalue true -command SetStrain ] + +#pack .figure .figure.scale + +# Create the canvas +set theCanvas [canvas .figure.canvas -bg "#ffffff" -height $height -width $width] + +$theCanvas create line 10 [expr $height/2.0] [expr $width-10] [expr $height/2.0] +$theCanvas create line [expr $width/2.0] 10 [expr $width/2.0] [expr $height-10] +$theCanvas create text [expr $width-30] [expr $height/2.0-10] -text "Strain $maxStrain" +$theCanvas create text [expr $width/2.0+30] 10 -text "Stress $maxStress" + +pack .figure .figure.canvas .figure.scale + +pack .figure -side top + +frame .style2 -borderwidth 1 -relief sunken +pack .style2 -fill x -pady 2 + +# ############################################################## +# Define the Settings frame & Settings procedure +# ############################################################## + +frame .settings -borderwidth 1 -relief raised + +label .settings.lstrain -text "Strain (max absolute)" +entry .settings.estrain -textvariable maxStrain -relief sunken +label .settings.lstress -text "Stress (max absolute)" +entry .settings.estress -textvariable maxStress -relief sunken +button .settings.ok -text OK -command "Reset" + +grid .settings.lstrain -row 0 -column 0 -sticky e +grid .settings.estrain -row 0 -column 1 -sticky ew +grid .settings.ok -row 0 -rowspan 2 -column 2 -sticky nsew +grid .settings.lstress -row 1 -column 0 -sticky e +grid .settings.estress -row 1 -column 1 -sticky ew + +# Procedure Settings - used to set bottom frame to be the settings frame +proc Settings { } { + global toggleFrame + global .settings + + pack forget $toggleFrame + set toggleFrame .settings + pack $toggleFrame -side bottom -fill x -pady 4 +} + + +# ############################################################## +# Define the Values frame & SetValues procedure +# ############################################################## + +frame .values -borderwidth 1 -relief raised + +label .values.strain -text "Strain : " +label .values.stress -text "Stress : " +label .values.tangent -text "Tangent: " +label .values.dstrain -text " 0.0" +label .values.dstress -text " 0.0" +label .values.dtangent -text " 0.0" + +grid .values.strain -row 0 -column 0 -sticky e +grid .values.dstrain -row 0 -column 1 -sticky ew +grid .values.stress -row 1 -column 0 -sticky e +grid .values.dstress -row 1 -column 1 -sticky ew +grid .values.tangent -row 2 -column 0 -sticky e +grid .values.dtangent -row 2 -column 1 -sticky ew + +# Procedure Settings - used to set bottom frame to be the settings frame +proc SetValues { } { + global toggleFrame + global .values + + pack forget $toggleFrame + set toggleFrame .values + pack $toggleFrame -side bottom -fill x -pady 4 +} + +# ############################################################## +# Define the Reset Procedure +# ############################################################## + +set pointTag 0 + +proc Reset { } { + global theCanvas + global theScale + global height + global width + global maxStrain + global maxStress + global strain + global xLast + global yLast + global width + global pointTag + + set strain 0 + + $theScale configure -from [expr -$maxStrain] -to $maxStrain -length $width \ + -tickinterval [expr $maxStrain / 2] -resolution [expr $maxStrain / 100] \ + -showvalue true -command SetStrain + + $theCanvas move all 400 400 + + $theCanvas create line 10 [expr $height/2.0] [expr $width-10] [expr $height/2.0] + $theCanvas create line [expr $width/2.0] 10 [expr $width/2.0] [expr $height-10] + $theCanvas create text [expr $width-30] [expr $height/2.0-10] -text "Strain $maxStrain" + $theCanvas create text [expr $width/2.0+30] 10 -text "Stress $maxStress" +# $theCanvas create rectangle 30 30 40 40 -tags "drawing" + + set pointTag [ + $theCanvas create oval [expr $height/2.-3] [expr $width/2.-3] [expr $height/2.+3] [expr $width/2.+3] \ + -tags "point" -fill "red" + ] +# $theCanvas create rectangle [expr $height/2.-2] [expr $width/2.-2] [expr $height/2.+2] [expr $width/2.+2] -tags "point" + + puts $pointTag + + set xLast [expr $width/2] + set yLast [expr $height/2] + + .values.dstrain config -text [format "%6.4e" 0.0] + .values.dstress config -text [format "%6.4e" 0.0] + .values.dtangent config -text [format "%6.4e" 0.0] + + # ResetMaterial +} + + +proc ResetMaterial { } { + global matID + if {$matID != 0} {eval uniaxialTest $matID} +} + + +# ############################################################## +# Define the SetStrain Procedure +# ############################################################## + +proc SetStrain {strain} { + global theCanvas + global maxStrain + global maxStress + global xLast + global yLast + global matID + global height + global width + global toggleFrame + global pointTag + + + if {$matID != 0} { + invoke UniaxialMaterial $matID "strain $strain" + invoke UniaxialMaterial $matID "commit" + set stress [invoke UniaxialMaterial $matID stress ] + set tangent [invoke UniaxialMaterial $matID tangent] + + puts "$strain $stress $tangent" + + set fileOpen [open results.out a] + puts $fileOpen "$strain $stress $tangent" + close $fileOpen + + set diffStrain [expr $width/(2*$maxStrain)] + set diffStress [expr $height/(2*$maxStress)] + + set x [expr $width / 2 + $strain * $diffStrain] + set y [expr $height / 2 - $stress * $diffStress] + $theCanvas create line $xLast $yLast $x $y + + $theCanvas move $pointTag [expr $x-$xLast] [expr $y-$yLast] + + set xLast $x + set yLast $y + + if {$toggleFrame == ".values"} { + .values.dstrain config -text [format "%6.4e" $strain] + .values.dstress config -text [format "%6.4e" $stress] + .values.dtangent config -text [format "%6.4e" $tangent] + } + } +} + +tkloop diff --git a/tools/SCRIPTS/toOpenSeesPy.py b/tools/SCRIPTS/toOpenSeesPy.py new file mode 100644 index 0000000000..75bbf542a7 --- /dev/null +++ b/tools/SCRIPTS/toOpenSeesPy.py @@ -0,0 +1,116 @@ +# Convert an OpenSees(Tcl) script to OpenSeesPy +# Author: Michael H. Scott, michael.scott@oregonstate.edu +# Date: June 2018 +# +# Usage in a Python script +# exec(open('toOpenSeesPy.py').read()) +# ... +# outfile = open('model.py','w') +# toOpenSeesPy('model.tcl',outfile) +# toOpenSeesPy('anotherScript.tcl',outfile) +# ... +# outfile.close() +# +# - Assumes the OpenSees(.tcl) file defines the model line by line +# without any loops, conditionals, variables, expressions, etc. +# This is the format generated when you export a model from +# OpenSees Navigator and perhaps from other front-ends to OpenSees. +# +# - The calling Python script should open and close the file stream for +# for writing the converted .py file. This allows you to call the +# converter on multiple Tcl files in sequence, as shown above. +# +# - If your OpenSees(.tcl) file uses any loops, conditionals, variables, +# expressions, etc., you might be better off to port your OpenSees +# model from Tcl to Python manually, or you can look in to Tkinter. +# +# - You may have some luck making your own "middleware" to convert your +# OpenSees(.tcl) script to a model defined line by line by inserting +# output statements in your loops and other constructs. Even though +# this won't get you to 100% conversion and you'll still have some +# conversions to make here and there, it'll get you pretty far. +# +# set output [open lineByLine.tcl w] +# ... +# for {set i 1} {$i <= $N} {incr i} { +# element truss $i ... +# puts $output "element truss $i ..." +# } +# ... +# close $output +# +# Then, in your Python script, call toOpenSeesPy with lineByLine.tcl as +# the input file. +# +# - If you see any improvements to make to this toOpenSeesPy function, +# please submit a pull request at OpenSees/OpenSees on github + + +# Helper function to deterine if a variable is a floating point value or not +# +def isfloat(value): + try: + float(value) + return True + except ValueError: + return False + +# Function that does the conversion +# +def toOpenSeesPy(infile, outfile): + outfile.write('\n\n') + infile = open(infile,'r') + for line in infile: + info = line.split() + N = len(info) + + # Ignore a close brace + if N > 0 and info[0][0] == '}': + continue + # Echo a comment line + if N < 2 or info[0][0] == '#': + outfile.write(line) + continue + + # Needs to be a special case for now due to beam integration + if info[1] == 'forceBeamColumn' or info[1] == 'dispBeamColumn': + secTag = info[6] + eleTag = info[2] + Np = info[5] + if info[1] == 'dispBeamColumn': + outfile.write('ops.beamIntegration(\'Legendre\',%s,%s,%s)\n' % (eleTag,secTag,Np)) + if info[1] == 'forceBeamColumn': + outfile.write('ops.beamIntegration(\'Lobatto\',%s,%s,%s)\n' % (eleTag,secTag,Np)) + outfile.write('ops.element(\'%s\',%s,%s,%s,%s,%s)\n' % (info[1],eleTag,info[3],info[4],info[7],eleTag)) + continue + + # Change print to printModel + if info[0] == 'print': + info[0] = 'printModel' + + # For everything else, have to do the first one before loop because of the commas + if isfloat(info[1]): + outfile.write('ops.%s(%s' % (info[0],info[1])) + else: + outfile.write('ops.%s(\'%s\'' % (info[0],info[1])) + # Now loop through the rest with preceding commas + writeClose = True + for i in range (2,N): + if info[i] == '{': + writeClose = True + break + if info[i] == '}': + writeClose = False + break + if info[0] == 'recorder': + # If it's a recorder, make everything immediately after material, section, or fiber a string + if info[i-1] in ['material','section','fiber'] and isfloat(info[i]): + outfile.write(',str(%s)' % info[i]) + continue + if isfloat(info[i]): + outfile.write(',%s' % info[i]) + else: + outfile.write(',\'%s\'' % info[i]) + if writeClose: + outfile.write(')\n') + infile.close() diff --git a/tools/SCRIPTS/uniaxialMaterials.tcl b/tools/SCRIPTS/uniaxialMaterials.tcl new file mode 100644 index 0000000000..3873e14c7c --- /dev/null +++ b/tools/SCRIPTS/uniaxialMaterials.tcl @@ -0,0 +1,27 @@ +# source for opstkUniaxialMaterialViewer.tcl +# define sample input, add to Material menu +set Es 29000 +set fy 50.0 +set ey [expr $fy/$Es] +set R0 15.0 +set b 0.05 +set cR1 0.925 +set cR2 0.150 +$m add command -label "Elastic" -command "uniaxialMaterialSample Elastic matTag=1 E=$Es" +$m add command -label "ElasticPP" -command "uniaxialMaterialSample ElasticPP matTag=1 E=$Es yieldStrain=$ey" +$m add command -label "Steel01" -command "uniaxialMaterialSample Steel01 matTag=1 Fy=$fy E=$Es b=$b" +$m add command -label "Steel02" -command "uniaxialMaterialSample Steel02 matTag=1 Fy=$fy E=$Es b=$b R0=$R0 CR1=$cR1 CR2=$cR2" +$m add command -label "Steel04" -command \ + "uniaxialMaterialSample Steel4 matTag=1 fy=$fy E=$Es -kin bk=$b R0=$R0 r1=$cR1 r2=$cR2"; # -iso bi=$bi $rho_i bl=$bl $R_i l_yp=1 < -ult $f_u $R_u > -mem 10 " +$m add command -label "BoucWenOriginal" -command "uniaxialMaterialSample BoucWenOriginal matTag=1 E=$Es fy=$fy alphaL=0.025 alphaNL=0.0 mu=2.0 eta=1.0 beta=0.55 gamma=0.45" +$m add command -label "BoucWen" -command "uniaxialMaterialSample BoucWen matTag=1 alpha=0.025 ko=41760.0 n=1.0 gamma=0.45 beta=0.55 Ao=1.0 deltaA=0.0 deltaNu=0.0 deltaEta=0.0" +$m add command -label "ReinforcingSteel" -command \ + "uniaxialMaterialSample ReinforcingSteel matTag=1 fy=$fy fu=120 Es=$Es Esh=20000 esh=0.1 eult=1.0" +$m add command -label "MultiLinear" -command "uniaxialMaterialSample MultiLinear matTag=1 e1=0.05 s1=10 e2=0.15 s2=20 e3=0.25 s3=25 e4=0.6 s4=0." +$m add command -label "SimpleFracture" -command "uniaxialMaterialSample SimpleFracture matTag=1 otherTag=1 maxStrain=0.4" +$m add command -label "MultiLinear2" -command "uniaxialMaterialSample MultiLinear matTag=1 e1=0.0001341737342935323 s1=146.400 e2=0.008252161648910277 s2=158.720000 e3=0.022357133345899175 s3=0.0 e4=10.0223571333459 s4=0." +$m add command -label "Bilin02" -command "uniaxialMaterialSample Bilin02 matTag=1 K0=3672850.0 asP=0.00082 asN=0.00082 myP=24640.0 myN=-24640.0 LS=1.60688 LC=1.60688 LA=1.60688 LK=1.60688 cS=1.0 cC=1.0 cA=1.0 cK=1.0 tpP=0.405899 tpN=0.405899 tpcP=0.705248 tpcN=0.705248 rP=0.4 rN=0.4 tuP=0.4 tuN=0.4 dP=1.0 dN=1.0 nFactor=1."; +$m add command -label "ModIMKPeak" -command \ + "uniaxialMaterialSample ModIMKPeakOriented matTag=1 K0=3672850.0 asP=0.00082 asN=0.00082 myP=24640.0 myN=-24640.0 LS=1.60688 LC=1.60688 LA=1.60688 LK=1.60688 cS=1.0 cC=1.0 cA=1.0 cK=1.0 tpP=0.405899 tpN=0.405899 tpcP=0.705248 tpcN=0.705248 rP=0.4 rN=0.4 tuP=0.4 tuN=0.4 dP=1.0 dN=1.0"; +$m add command -label "Bilin" -command "uniaxialMaterialSample Bilin matTag=1 p1=3672850.0 p2=0.0008263977213611486 p3=0.0008263977213611486 p4=24640.0 p5=-24640.0 p6=1.60688 p7=1.60688 p8=1.6068826060003187 p9=1.60688 p10=1.0 p11=1.0 p12=1.0 p13=1.0 p14=0.405899 p15=0.405899 p16=0.705248 p17=0.705248 p18=0.4 p19=0.4 p20=0.4 p21=0.4 p22=1.0 p23=1.0 p24=1."; +$m add command -label "HoehlerStanton" -command "uniaxialMaterialSample HoehlerStanton matTag=1 200000 0.015 400 550 0.15 3 25 23 0.2" diff --git a/tools/apidoc.lua b/tools/apidoc.lua new file mode 100644 index 0000000000..a046b7f849 --- /dev/null +++ b/tools/apidoc.lua @@ -0,0 +1,30 @@ + +-- pandoc's List type +local List = require 'pandoc.List' + + +--- Filter function for code blocks +function Div (cb) + local apidoc = cb.attributes['apidoc'] + if not apidoc then + return + end + io.stderr:write("apidoc.lua: " .. apidoc .. "\n") + + local blocks = List:new() + fh = io.popen("python -m opensees.emit.apidoc " .. apidoc) + local contents = pandoc.read(fh:read '*all', "markdown").blocks + fh:close() + -- contents = pandoc.walk_block( + -- pandoc.Div(contents), { } + -- ).content + blocks:extend(contents) + -- blocks:extend(cb) + -- return blocks + return blocks +end + + + + + diff --git a/tools/build-conda-package.sh b/tools/build-conda-package.sh new file mode 100644 index 0000000000..c23028e03b --- /dev/null +++ b/tools/build-conda-package.sh @@ -0,0 +1,52 @@ +#!/bin/bash + +# SOURCE: +# https://giswqs.medium.com/building-a-conda-package-and-uploading-it-to-anaconda-cloud-6a3abd1c5c52 + +# change the package name to the existing PyPi package you would like to build +# and adjust the Python versions + +pkg='opensees' +recipe="etc/conda" +python_versions=( 3.7 3.8 3.9 ) +platforms=( osx-64 linux-32 linux-64 win-32 win-64 ) + +alias CONDA-BUILD="conda mambabuild" + +if "" then + echo "Building conda package ..." + cd ~ + conda skeleton pypi $pkg + cd $pkg + wget https://conda.io/docs/_downloads/build1.sh + wget https://conda.io/docs/_downloads/bld.bat + cd ~ +fi + +# building conda packages +for v in "${python_versions[@]}" +do + CONDA-BUILD --python $v $recipe +done + +# convert package to other platforms +cd ~ +find $HOME/*/conda-bld/linux-64/ -name *.tar.bz2 | while read file +do + echo "converting $file" + #conda convert --platform all $file -o $HOME/conda-bld/ + for platform in "${platforms[@]}" + do + conda convert --platform $platform $file -o $HOME/conda-bld/ + done + +done + +# upload packages to conda +find $HOME/conda-bld/ -name *.tar.bz2 | while read file +do + echo "uploading $file" + anaconda upload $file +done + +echo "Building conda package complete!" diff --git a/tools/doc.py b/tools/doc.py new file mode 100644 index 0000000000..14a43a1aca --- /dev/null +++ b/tools/doc.py @@ -0,0 +1,173 @@ +import re +import sys +import inspect +import textwrap +import importlib + + +from opensees import ast, obj + + +write = lambda *args, **kwds: print(*args, **kwds) + + +def get_type(a): + pass + +def write_grp(a): + args = (write_grp(i) if isinstance(i, ast.Grp) else (i.name if i.name else "") for i in a.args) + return "["+",".join(args)+"]" + #write(f"`{name}``typ`{about}", end="") + +def write_arg(a): + write("",end="") + typ = f"{a.__class__.__name__}" + # about = a.about.replace('\n',' ').strip() + name = a.name or "" + default = " = "+str(a.default) if a.default else "" + about = re.sub('[\s+]', ' ', a.about.replace('\n',' ')) + # write(f">| {a.name} | `{typ}` | {about} |") + # print(f">>>> {name} ") + + if isinstance(a, ast.Grp): + if name == "": + name = f"{write_grp(a)}" + #elif a.type: + # typ = f"{a.type.__name__}" + else: + typ = f"{write_grp(a)}" + write(f"{name+default}{typ}{about}", end="") + write("") + [write_arg(i) for i in a.args] + write("
") + + elif isinstance(a, ast.Map): + # typ = f"{a.__class__.__name__}({a.key_type.__class__.__name__}: {a.val_type.__name__})" + typ = f"{a.__name__}" + write(f"{name+default}{typ}{about}", end="") + write("") + write_arg(a.key_type) + write_arg(a.val_type) + write("
") + + else: + if isinstance(a, ast.Ref): + typ = f"{a.__class__.__name__}({a.type.__name__})" + write(f"{name+default}{typ}{about}", end="") + + write("") + +def write_obj(v): + s = str(inspect.signature(v)).replace('=None','') + if len(s) > 65: + # + s = s.replace(", ", ",
   ") + write(textwrap.dedent(f""" + + + {k}{s} + + """)) + try: + write(v.about) + except: + write(v.__doc__ or "") + if hasattr(v, "_img"): + write(f"![](/figures/{v._img})") + + # write(">| | | |\n|--|--|--|") + # write("") + write(textwrap.dedent(""" +
+ + + + + """)) + # write("") + try: + pass + except Exception as e: + print(e, file=sys.stderr) + if True: + for a in v._args: + # print(a, file=sys.stderr) + write_arg(a) + + write(textwrap.dedent(""" + +
+ + """)) + #write(":::\n") + + + +def write_single(obj_name): + lib = vars(getattr(lib, obj_name.split(".")[0])) + write(textwrap.dedent(getattr(lib, attrs[0]).__doc__ or "")) + + +if __name__=="__main__": + import opensees.lib + + + libs = [] + attrs = [] + solo = False + args = iter(sys.argv[1:]) + for arg in args: + if arg[0] != "-": + libs.append(arg) + elif arg == "--attr": + attrs.append(next(args)) + elif arg =="--single": + solo = True + + if solo: + lib = importlib.import_module(libs[0]) + write_single(lib, attrs[0]) + sys.exit() + + for lib_name in libs: + lib = importlib.import_module(lib_name) + head = (attrs or lib_name.split("."))[-1] + write(textwrap.dedent(f""" + --- + title: {head} + ... + + + + # {head} + + """)) + if not attrs: + objs = vars(lib) + write(textwrap.dedent(lib.__doc__ or "")) + + else: + objs = vars(getattr(lib, attrs[0])) + write(textwrap.dedent(getattr(lib, attrs[0]).__doc__ or "")) + + write('
') + + for k,v in objs.items(): + if isinstance(v, type): + if issubclass(v, obj.Component) and hasattr(v, "_args"): + write_obj(v) + else: + pass + write("
") + + + + + + diff --git a/tools/docopt_c.py b/tools/docopt_c.py new file mode 100644 index 0000000000..906705196f --- /dev/null +++ b/tools/docopt_c.py @@ -0,0 +1,701 @@ +#!/usr/bin/env python +# -*- coding:utf-8 -*- + +# Copyright (c) 2012 Vladimir Keleshev, +# (see LICENSE-MIT file for copying) + + +"""Usage: docopt_c.py [options] [] + +Processes a docopt formatted string, from either stdin or a file, and +outputs the equivalent C code to parse a CLI, to either the stdout or a file. + +Options: + -o, --output-name= + Filename used to write the produced C file. + If not present, the produced code is printed to stdout. + -t, --template=