From ca47334e9f918fd2919c3a8a9c50f0cdab79ec9a Mon Sep 17 00:00:00 2001 From: claudio perez <50180406+claudioperez@users.noreply.github.com> Date: Tue, 20 Aug 2024 15:36:33 -0700 Subject: [PATCH 1/8] cmp - add `surface` command --- src/opensees/openseespy.py | 665 ++++++++++++++++++++++--------------- 1 file changed, 404 insertions(+), 261 deletions(-) diff --git a/src/opensees/openseespy.py b/src/opensees/openseespy.py index 5eb657db6e..ba443eafcc 100644 --- a/src/opensees/openseespy.py +++ b/src/opensees/openseespy.py @@ -13,257 +13,11 @@ import json from functools import partial -from .tcl import Interpreter +from .tcl import Interpreter, _lift # something to compare the output of model.analyze to: successful = 0 -# 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", - "mesh" -} - class OpenSeesError(Exception): pass @@ -316,21 +70,62 @@ def _split_iter(source, sep=None, regex=False): 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) -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 + 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") @@ -389,6 +184,29 @@ def eval(self, cmd: str) -> str: 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): @@ -487,6 +305,7 @@ 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 @@ -522,6 +341,8 @@ def _mesh_line(self, tag, numnodes, ndtags, id, ndf:int, meshsize, eleType='', e if i < nn: 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): @@ -566,6 +387,9 @@ def __init__(self, *args, echo_file=None, **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() @@ -573,6 +397,15 @@ def asdict(self): def setFactor(self, factor): pass + def element(self, type, tag, *args, **kwds): + if tag is None: + tag = 0 + for existing_tag in self.getEleTags(): + 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") @@ -584,12 +417,78 @@ def getTangent(self, **kwds): 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 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 @@ -613,14 +512,258 @@ def _as_str_arg(arg, name: str = None): return str(arg) -class FedeasModel(Model): - @property - def nf(self): - return self.numDOF() # 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: From 667e099acb5d7fc60ae1747d11525b5b8c46f06f Mon Sep 17 00:00:00 2001 From: claudio perez <50180406+claudioperez@users.noreply.github.com> Date: Tue, 20 Aug 2024 15:41:06 -0700 Subject: [PATCH 2/8] cmp - error reporting --- src/opensees/__main__.py | 11 ++++++++++- src/opensees/tcl.py | 7 +++++-- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/opensees/__main__.py b/src/opensees/__main__.py index c343dc394c..045662ad98 100755 --- a/src/opensees/__main__.py +++ b/src/opensees/__main__.py @@ -1,4 +1,11 @@ #!/usr/bin/env python3 +""" +This file implements the primary command line interface for +the package which is invoked by running: + + python -m opensees + +""" import os import sys @@ -163,13 +170,15 @@ def parse_args(args): for cmd in opts["commands"]: tcl.eval(cmd) + # Store the exit code of the Tcl script code = 0 + # if script is not None: try: code = tcl.eval(script) - except opensees.tcl.tkinter._tkinter.TclError as e: + except opensees.tcl.TclError as e: print(e, file=sys.stderr) if not opts["interact"]: sys.exit(-1) diff --git a/src/opensees/tcl.py b/src/opensees/tcl.py index ec75e95103..4f1617d962 100644 --- a/src/opensees/tcl.py +++ b/src/opensees/tcl.py @@ -15,6 +15,8 @@ from opensees.library.obj import Component +class TclError(Exception): + pass def exec(script: str, silent=False, analysis=True)->dict: """ @@ -146,8 +148,9 @@ def eval(self, string): return self._tcl.tk.eval(string) except tkinter._tkinter.TclError as e: - print(self._tcl.getvar("errorInfo"), file=sys.stderr) - raise e + raise TclError(self._tcl.getvar("errorInfo")) +# print(self._tcl.getvar("errorInfo"), file=sys.stderr) +# raise e def serialize(self)->dict: import tempfile From f83da8489a599943ef70f5b5f529bc4194670bda Mon Sep 17 00:00:00 2001 From: claudio perez <50180406+claudioperez@users.noreply.github.com> Date: Tue, 20 Aug 2024 15:42:02 -0700 Subject: [PATCH 3/8] cmp - parallelize torsion solver --- src/opensees/section/torsion.py | 121 ++++++++++++++++++++++++-------- 1 file changed, 91 insertions(+), 30 deletions(-) diff --git a/src/opensees/section/torsion.py b/src/opensees/section/torsion.py index 0c74a2c005..4432fd932b 100644 --- a/src/opensees/section/torsion.py +++ b/src/opensees/section/torsion.py @@ -9,39 +9,78 @@ # https://civil-terje.sites.olt.ubc.ca/files/2020/02/Screenshot-Solid-Cross-section.pdf # import numpy as np +from functools import partial class TorsionAnalysis: pass -def assembleStiffnessMatrix(Ka, Kg, nodeList, nde, ndf): - nne = len(nodeList) +def assemble_matrix(Ka, ke, conn, nde, ndf): + nne = len(conn) for j in range(nne): for k in range(j + 1): for m in range(nde): for n in range(nde): - Ka[(nodeList[j] - 1)*ndf + m, (nodeList[k] - 1)*ndf + n] += Kg[j*nde + m, k*nde + n] + Ka[(conn[j] - 1)*ndf + m, (conn[k] - 1)*ndf + n] += ke[j*nde + m, k*nde + n] if j != k: - Ka[(nodeList[k] - 1)*ndf + n, (nodeList[j] - 1)*ndf + m] += Kg[j*nde + m, k*nde + n] - - - # This could have been done, less efficiently, by using Tga, like this: - # Tga = np.zeros((numDOFsInElement, nf)) - # for j in range(nodeList.size): - # for k in range(nde): - # Tga[j * nde + k, (nodeList[j] - 1) * ndf + k] = 1 - # Ka += (np.transpose(Tga).dot(Kg)).dot(Tga) + Ka[(conn[k] - 1)*ndf + n, (conn[j] - 1)*ndf + m] += ke[j*nde + m, k*nde + n] return Ka -def assembleLoadVector(Fa, Fg, nodes, nde, ndf): - nne = len(nodes) - for j in range(nne): +def assemble_vector(Fa, fe, nodes, nde, ndf): + nen = len(nodes) + for j in range(nen): for m in range(nde): - Fa[(nodes[j]-1) * ndf + m] += Fg[j * nde + m] + Fa[(nodes[j]-1) * ndf + m] += fe[j * nde + m] return Fa + +class TorsionTriangle: + _class_mass = 1/12*(np.eye(3) + np.ones((3,3))) + + def __init__(self, xyz): + ((y1, y2, y3), (z1, z2, z3)) = xyz + self.xyz = xyz + self.area = area = 0.5 * ((y2 - y1) * (z3 - z1) - (y3 - y1) * (z2 - z1)) + + z12 = z1 - z2 + z23 = z2 - z3 + z31 = z3 - z1 + y32 = y3 - y2 + y13 = y1 - y3 + y21 = y2 - y1 + + area = 0.5 * ((y2 - y1) * (z3 - z1) - (y3 - y1) * (z2 - z1)) + + k11 = ( y32**2 + z23**2) + k12 = (y13*y32 + z23*z31) + k13 = (y21*y32 + z12*z23) + k22 = ( y13**2 + z31**2) + k23 = (y13*y21 + z12*z31) + k33 = ( y21**2 + z12**2) + ke = 1/(4.0*area)*np.array([[k11, k12, k13], + [k12, k22, k23], + [k13, k23, k33]]) + + fe = -1/6.*np.array([ + ((y1*y32 - z1*z23) + (y2*y32 - z2*z23) + (y3*y32 - z3*z23)), + ((y1*y13 - z1*z31) + (y2*y13 - z2*z31) + (y3*y13 - z3*z31)), + ((y1*y21 - z1*z12) + (y2*y21 - z2*z12) + (y3*y21 - z3*z12))]) + + self.load = fe + self.stiffness = ke + + def getMass(self): + rho = 1.0 + return rho*self.area*self._class_mass + + def getStiffness(self): + return self.stiffness + + def getLoad(self): + return self.load + def torsion_element(xyz): ((y1, y2, y3), (z1, z2, z3)) = xyz @@ -60,6 +99,9 @@ def torsion_element(xyz): k22 = ( y13**2 + z31**2) k23 = (y13*y21 + z12*z31) k33 = ( y21**2 + z12**2) + + me = area/12*(np.eye(3) + np.ones((3,3))) + ke = 1/(4.0*area)*np.array([[k11, k12, k13], [k12, k22, k23], [k13, k23, k33]]) @@ -69,23 +111,42 @@ def torsion_element(xyz): ((y1*y13 - z1*z31) + (y2*y13 - z2*z31) + (y3*y13 - z3*z31)), ((y1*y21 - z1*z12) + (y2*y21 - z2*z12) + (y3*y21 - z3*z12))]) - return ke, fe + return me, ke, fe + +def _wrap_elem(centroid, nodes, conn): + return conn, torsion_element((nodes[conn] - centroid).T) +import multiprocessing def solve_torsion(section, mesh): ndf = 1 nde = 1 nodes = mesh.points nf = ndf*len(nodes) Ka = np.zeros((nf, nf)) + Ma = np.zeros((nf, nf)) Fa = np.zeros(nf) connectivity = mesh.cells[1].data - for conn in connectivity: - ((y1, y2, y3), (z1, z2, z3)) = (nodes[conn] - section.centroid).T - ke, fe = torsion_element((nodes[conn] - section.centroid).T) - nodeList = conn + 1 - Ka = assembleStiffnessMatrix(Ka, ke, nodeList, nde, ndf) - Fa = assembleLoadVector(Fa, fe, nodeList, nde, ndf) + + threads = 6 + chunk = 200 + with multiprocessing.Pool(threads) as pool: + for conn, (me, ke, fe) in pool.imap_unordered( + partial(_wrap_elem, section.centroid, nodes), + mesh.cells[1].data, + chunk + ): + Ka = assemble_matrix(Ka, ke, conn+1, nde, ndf) + Ma = assemble_matrix(Ma, me, conn+1, nde, ndf) + Fa = assemble_vector(Fa, fe, conn+1, nde, ndf) + +# queue = ((conn, torsion_element((nodes[conn] - section.centroid).T)) +# for conn in connectivity) + +# for conn, (me, ke, fe) in queue: +# Ka = assemble_matrix(Ka, ke, conn+1, nde, ndf) +# Ma = assemble_matrix(Ma, me, conn+1, nde, ndf) +# Fa = assemble_vector(Fa, fe, conn+1, nde, ndf) # Lock the solution at one node and solve for the others Pf = Fa[:nf-1] @@ -98,8 +159,8 @@ def solve_torsion(section, mesh): #return ua -#def shear_centre(section, mesh, ua): - # Determine shear centre coordinates +#def shear_center(section, mesh, ua): + # Determine shear center coordinates A = section.area Iy = section.iyc Iz = section.ixc @@ -118,7 +179,7 @@ def solve_torsion(section, mesh): #def normalize_shear(): thetaPrincipal = 0.0 - # Normalizing constant and shear centre coordinates + # Normalizing constant and shear center coordinates # normalizingConstant = -warpIntegral / A if np.abs(thetaPrincipal) > 1e-3: ysc = -yscIntegral / IyRotated @@ -243,11 +304,11 @@ def torsion_constant(section, mesh, warp, sc): # print('\n'"Product of inertia Iyz: %.2f" % Iyz) # print('\n'"New principal moment of inertia Iy: %.2f" % IyRotated) # print('\n'"New principal moment of inertia Iz: %.2f" % IzRotated) -# print('\n'"Shear centre y-coordinate in original coordinate system: %.2f" % ysc) -# print('\n'"Shear centre z-coordinate in original coordinate system: %.2f" % zsc) +# print('\n'"Shear center y-coordinate in original coordinate system: %.2f" % ysc) +# print('\n'"Shear center z-coordinate in original coordinate system: %.2f" % zsc) # else: -# print('\n'"Shear centre y-coordinate: %.2f" % ysc) -# print('\n'"Shear centre z-coordinate: %.2f" % zsc) +# print('\n'"Shear center y-coordinate: %.2f" % ysc) +# print('\n'"Shear center z-coordinate: %.2f" % zsc) # print('\n'"St. Venant torsion constant J: %.2f" % J) # print('\n'"Warping torsion constant Cw: %.2f" % Cw) # From c3998d3ba3441cab7bdd30ed4abe7da3d7b97267 Mon Sep 17 00:00:00 2001 From: claudio perez <50180406+claudioperez@users.noreply.github.com> Date: Tue, 20 Aug 2024 15:42:37 -0700 Subject: [PATCH 4/8] cmp - organize section utilities --- src/opensees/section/patch.py | 2 +- src/opensees/section/section.py | 349 ++++---------------------------- 2 files changed, 37 insertions(+), 314 deletions(-) diff --git a/src/opensees/section/patch.py b/src/opensees/section/patch.py index f5b64738de..27f0dd5464 100644 --- a/src/opensees/section/patch.py +++ b/src/opensees/section/patch.py @@ -201,7 +201,7 @@ def fibers(self): @_patch class quad(_Polygon): """A quadrilateral shaped patch. - The geometry of the patch is defined by four vertices: I J K L. + The geometry of the patch is defined by four vertices: I J K L. The coordinates of each of the four vertices is specified in *counter clockwise* sequence """ _img = "quadPatch.svg" diff --git a/src/opensees/section/section.py b/src/opensees/section/section.py index b2c98db1cb..da98b2e3be 100644 --- a/src/opensees/section/section.py +++ b/src/opensees/section/section.py @@ -7,27 +7,45 @@ from math import pi, sin, cos from opensees.library import LibCmd, Cmd, Component from opensees.library import uniaxial +from opensees.library.frame import Area, Iyc, Ixc, Yng from opensees.library.ast import Tag, Num, Blk, Ref, Flg, Map from . import patch -layer = patch.layer _section = LibCmd("section") class _FiberCollection: - def __init__(self, areas): - self.areas = areas + def __init__(self, shapes): + self.shapes = shapes def __contains__(self, point): - return any(point in area for area in self.areas) + return any(point in area for area in self.shapes) + +@_section +class ElasticSection: +# section Elastic $secTag $E $A $Iz <$G $alphaY> +# section Elastic $secTag $E $A $Iz $Iy $G $J <$alphaY $alphaZ> + _ndms = (2, 3) + _args = [ + Tag(), + Yng(), + Area(), + Ixc(), + Iyc(), + Num("G", about="Shear Modulus (optional for 2D analysis, required for 3D analysis)", ndm_reqd=(3,), reqd=False), + Num("J", about="torsional moment of inertia of section (required for 3D analysis)", ndm_reqd=(3,), reqd=False), + Num("alphaY", about="shear shape factor along the local $y$-axis", reqd=False), + Num("alphaZ", about="shear shape factor along the local $z$-axis", reqd=False), + ] @_section class FiberSection(_FiberCollection): + _ndms = (2, 3) _args = [ Tag(), - Num("GJ", flag="-GJ", field="torsional_stiffness", optional=True, - about="linear-elastic torsional stiffness assigned to the section (optional, default = no torsional stiffness)"), - Blk("areas", from_prop="fibers", default=[], type=Cmd, defn=dict( + Num("GJ", flag="-GJ", field="torsional_stiffness", ndm_reqd=(3,), optional=True, + about="linear-elastic torsional stiffness $GJ$"), + Blk("shapes", from_prop="fibers", default=[], type=Cmd, defn=dict( fiber=LibCmd("fiber"), layer=LibCmd("layer") ) @@ -43,7 +61,7 @@ def init(self): self._area = None if "material" in self.kwds: mat = self.kwds["material"] - for i in self.areas: + for i in self.shapes: if i.material is None: i.material = mat @@ -56,11 +74,11 @@ def get_refs(self): def add_patch(self, patch): self._fibers = None - self.areas.append(patch) + self.shapes.append(patch) def add_patches(self, patch): self._fibers = None - self.areas.extend(patch) + self.shapes.extend(patch) def __repr__(self): import textwrap @@ -73,17 +91,17 @@ def __repr__(self): @property def patches(self): - return [p for p in self.areas if p.get_cmd()[0] == "patch"] + return [p for p in self.shapes if p.get_cmd()[0] == "patch"] @property def layers(self): - return [p for p in self.areas if p.get_cmd()[0] == "layer"] + return [p for p in self.shapes if p.get_cmd()[0] == "layer"] @property def fibers(self): if self._fibers is None: self._fibers = [ - f for a in (a.fibers if hasattr(a,"fibers") else [a] for a in self.areas) + f for a in (a.fibers if hasattr(a,"fibers") else [a] for a in self.shapes) for f in a ] return self._fibers @@ -91,20 +109,20 @@ def fibers(self): @property def area(self): if self._area is None: - self._area = sum(i.area for i in self.areas) + self._area = sum(i.area for i in self.shapes) return self._area @property def centroid(self): # TODO: cache - return sum(i.centroid * i.area for i in self.areas) / self.area + return sum(i.centroid * i.area for i in self.shapes) / self.area @property def ixc(self): # TODO: cache yc = self.centroid[1] return sum( - p.ixc + (p.centroid[1]-yc)**2*p.area for p in self.areas + p.ixc + (p.centroid[1]-yc)**2*p.area for p in self.shapes ) @property @@ -112,7 +130,7 @@ def iyc(self): # TODO: cache xc = self.centroid[0] return sum( - p.iyc + (p.centroid[0]-xc)**2*p.area for p in self.areas + p.iyc + (p.centroid[0]-xc)**2*p.area for p in self.shapes ) @@ -121,7 +139,7 @@ def moic(self): # TODO: cache return [ [p.moi[i] + p.centroid[i]**2*p.area for i in range(2)] + [p.moi[-1]] - for p in self.areas + for p in self.shapes ] def print_properties(self): @@ -131,151 +149,6 @@ def print_properties(self): Iyy (centroid) | {self.iyc:9.0f} """)) - -def PolygonRing(n, extRad, intRad): - """ - Create a polygon annulus. - """ - psi = 2*pi/n - phi = psi - collection = [] - cover_divs = 1,2 # divisions in each slice of the cover - iR1, iR2 = [intRad/cos(pi/n)]*2 - oR1, oR2 = [extRad/cos(pi/n)]*2 - j = 0 - for i in range(n): - startAngle = (i - 1/2)*psi - sita1 = startAngle + j*phi # Slice start angle - sita2 = sita1 + phi # Slice end angle - # Cover Patch connects the circular core to the polygonal cover - collection.append( - patch.quad(None, cover_divs, - vertices = [ - [ iR1*cos(sita1), iR1*sin(sita1)], - [ oR1*cos(sita1), oR1*sin(sita1)], - [ oR2*cos(sita2), oR2*sin(sita2)], - [ iR2*cos(sita2), iR2*sin(sita2)], - ] - )) - sect = FiberSection(areas=collection) - sect.extRad = extRad - sect.intRad = intRad - return sect - - -def _oct_outline(Rcol): - n = 8 - phi = 2*pi/n - R = Rcol/cos(phi/2) - region = [ - layer.line(vertices=[ - [R*cos(i*phi-phi/2), R*sin(i*phi-phi/2)], - [R*cos(i*phi+phi/2), R*sin(i*phi+phi/2)] - ]) for i in range(n) - ] - sect = FiberSection(areas=region) - sect.extRad = Rcol - sect.intRad = 0.0 - return sect - -def ConfiningPolygon(n, extRad=None, intRad=None, divs=None, diameter=None, s=4, material=None): - psi = 2*pi/n - phi = psi/s - collection = [] - if divs is None: divs = 8,8 # divisions in each slice of the cover - iR1, iR2 = [intRad]*2 - - for i in range(n): - startAngle = (i - 1/2)*psi - for j in range(s): - sita1 = startAngle + j*phi # Slice start angle - sita2 = sita1 + phi # Slice end angle - oR1 = extRad/cos(pi/n - j*phi) - oR2 = extRad/cos(pi/n - (j+1)*phi) - # Cover Patch connects the circular core to the polygonal cover - collection.append( - patch.quad(None, divs, - vertices = [ - [ iR1*cos(sita1), iR1*sin(sita1)], - [ oR1*cos(sita1), oR1*sin(sita1)], - [ oR2*cos(sita2), oR2*sin(sita2)], - [ iR2*cos(sita2), iR2*sin(sita2)], - ] - )) - sect = FiberSection(areas=collection, material=material) - sect.extRad = extRad - sect.intRad = intRad - return sect - -def ConfinedPolygon( - n: int, - extRad: float, - intRad: float = None, - sdivs : tuple = None, # sector divs - s: int = 4, # number of patches per edge - DLbar = 4, - core_conc = None, - cover_conc = None, - ColMatTag = None, - units = None, - diameter = None -): - """ - Dcol : Width of octagonal column (to flat sides) - nLbar : Number of longitudinal bars - DLbar : Diameter of longitudinal bars - sTbar : Spacing of transverse spiral reinforcement - """ - # - # Column component dimensions - # - if intRad == extRad: - return _oct_outline(extRad) - - if intRad is None: - intRad = extRad - 2.0 - - # assert intRad == 0.0 and extRad > 0.0 - - inch = 1.0 - Dcol = 2*extRad - # tcover = 2.0*inch # 2 inch cover width - Rcol = Dcol/2.0 # Radius of octagonal column (to flat sides) - Dcore = 2*intRad - Rcore = Dcore/2.0 # Radius of circular core - # Along = pi*DLbar**2/4.0 # Area of longitudinal reinforcement bar - DTbar = 0.625*inch # Diameter of transverse spiral bar (#5 Rebar) - Asp = pi*DTbar**2/4.0 # Area of transverse spiral reinforcement bar - Dtran = Dcore - DTbar # Diameter of spiral of transverse spiral reinforcement - - # Density of transverse spiral reinforcement - # rho = 4.0* Asp/(Dtran*sTbar) - # Diameter of ring of longitudinal reinforcement - Dlong = Dcore - 2*DTbar - DLbar - # Rlong = Dlong/2. # Radius of ring of longitudinal reinforcement - - # Build Octagonal RC Column Section - cdivs = 64 # # fibers around the entire circumference - - numSlices = 8 # # slices in each of the 8 sections of the polygon - #if sdivs is None: sdivs = [1, 2] - sect = FiberSection( - material = ColMatTag, - GJ = 1e12, - areas = [ - # Inner Core Patch, 5 radial fibers - patch.circ(core_conc, [cdivs, 5], [0., 0.], 0.0, Rcore/2, 0.0, 2*pi), - # Outer Core Patch, 10 radial fibers - patch.circ(core_conc, [cdivs, 10], [0., 0.], Rcore/2, Rcore, 0.0, 2*pi) - ]) - - sect.add_patches(ConfiningPolygon(n, extRad, intRad, divs=sdivs, s=s).patches) - - #layer circ long_steel nLbar Along 0. 0. Rlong; # Longitudinal Bars - sect.extRad = extRad - sect.intRad = intRad - return sect - @_section class SectionAggregator: """ @@ -321,153 +194,3 @@ class SectionAggregator: authors=["Micheal H. Scott"] - - -def sect2shapely(section): - """ - Generate `shapely` geometry objects - from `opensees` patches or a FiberSection. - """ - import numpy as np - import shapely.geometry - from shapely.ops import unary_union - shapes = [] - if hasattr(section, "patches"): - patches = section.patches - else: - patches = [section] - - for patch in patches: - name = patch.__class__.__name__.lower() - if name in ["quad", "poly", "rect", "_polygon"]: - points = np.array(patch.vertices) - width,_ = points[1] - points[0] - _,height = points[2] - points[0] - shapes.append(shapely.geometry.Polygon(points)) - else: - n = 64 - x_off, y_off = 0.0, 0.0 - # calculate location of the point - external = [[ - 0.5 * patch.extRad * np.cos(i*2*np.pi*1./n - np.pi/8) + x_off, - 0.5 * patch.extRad * np.sin(i*2*np.pi*1./n - np.pi/8) + y_off - ] for i in range(n) - ] - if patch.intRad > 0.0: - internal = [[ - 0.5 * patch.intRad * np.cos(i*2*np.pi*1./n - np.pi/8) + x_off, - 0.5 * patch.intRad * np.sin(i*2*np.pi*1./n - np.pi/8) + y_off - ] for i in range(n) - ] - shapes.append(shapely.geometry.Polygon(external, [internal])) - else: - shapes.append(shapely.geometry.Polygon(external)) - - if len(shapes) > 1: - return unary_union(shapes) - else: - return shapes[0] - -def sect2gmsh(sect, size, **kwds): - import pygmsh - import numpy as np - if isinstance(size, int): - size = [size]*2 - - shape = sect2shapely(sect) - - with pygmsh.geo.Geometry() as geom: - geom.characteristic_length_min = size[0] - geom.characteristic_length_max = size[1] - coords = np.array(shape.exterior.coords) - holes = [ - geom.add_polygon(np.array(h.coords)[:-1], size[0], make_surface=False).curve_loop - for h in shape.interiors - ] - if len(holes) == 0: - holes = None - - poly = geom.add_polygon(coords[:-1], size[1], holes=holes) - # geom.set_recombined_surfaces([poly.surface]) - mesh = geom.generate_mesh(**kwds) - - mesh.points = mesh.points[:,:2] - for blk in mesh.cells: - blk.data = blk.data.astype(int) - # for cell in mesh.cells: - # cell.data = np.roll(np.flip(cell.data, axis=1),3,1) - return mesh - - -class TorsionalConstantAnalysis: - pass - -class MomentCurvatureAnalysis: - @staticmethod - def solve_eps(sect, kap, axial: float, eps0, tol=1e-6, maxiter=25): - # Newton-Raphson iteration - eps = eps0 - s = sect.getStressResultant([eps, kap], False) - for i in range(maxiter): - if abs(s[0] - axial) < tol: - return eps - s = sect.getStressResultant([eps, kap], False) - eps -= (s[0] - axial)/sect.getSectionTangent()[0,0] - return eps - - def __init__(self, axial): - pass - -def MomentSearch(section, ): - pass - -class MomentAxialLocus: - def __init__(self, section, axial): - self.axial = axial - self.section = section - - def plot(self): - pass - - def analyze(self, nstep = 30, incr=5e-6): - import matplotlib.pyplot as plt - import numpy as np - fig, ax = plt.subplots(1,2, constrained_layout=True) - sect = self.section - axial = self.axial - - if sect.name is None: - sect.name = 1 - - solve_eps = MomentCurvatureAnalysis.solve_eps - - # Curvature increment - dkap = incr - for P in axial: - with sect as s: - k0 = 0.0 - e0 = solve_eps(s, k0, P, solve_eps(s, k0, P, 0.0)) - PM = [ - s.getStressResultant([e0, k0], True), - s.getStressResultant([solve_eps(s, k0+dkap, P, e0), k0+dkap], True), - ] - e = e0 - kap = 2*dkap - for _ in range(nstep): - if abs(PM[-1][1]) < 0.995*abs(PM[-2][1]): - break - e = solve_eps(s, kap, P, e) - PM.append(s.getStressResultant([e, kap], True)) - kap += dkap - - p, m = zip(*PM) - - ax[0].plot(np.linspace(0.0, kap, len(m)), m) - - ax[1].scatter(m, p, s=0.2, color="k") - - ax[1].set_ylabel("Axial force, $P$") - ax[1].set_xlabel("Moment, $M$") - - plt.show() - From 15c1e59717649eae5c81f4309bf6f7cc814c6632 Mon Sep 17 00:00:00 2001 From: claudio perez <50180406+claudioperez@users.noreply.github.com> Date: Tue, 20 Aug 2024 15:43:02 -0700 Subject: [PATCH 5/8] cmp - clean up, organize --- src/opensees/library/model.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/opensees/library/model.py b/src/opensees/library/model.py index 3639592c17..22e8878b8f 100644 --- a/src/opensees/library/model.py +++ b/src/opensees/library/model.py @@ -103,8 +103,6 @@ def block(self, divs, type: str, args, points, **kwds): Ref(f"Node{1+i}", type=Node, attr="name", about="") for i in range(4) ]) - arg_spec = [Tag(), nodes]+[Str(f"a{i+1}") for i in range(len(args))] - elem_class = Ele(type, type, args=arg_spec) if len(self.m_nodes) > 0: join = dict(nodes={int(v.name): v.crd for k,v in self.m_nodes.items()}, cells=self.m_elems.keys()) @@ -113,20 +111,23 @@ def block(self, divs, type: str, args, points, **kwds): join = None - nodes, elems = shps.block.block(divs, block_type, points=points, + nodes, cells = shps.block.block(divs, block_type, points=points, append=False, join=join, **kwds) - for tag, coord in nodes.items(): self.node(tag, *coord) - for tag, nodes in elems.items(): + arg_spec = [Tag(), nodes]+[Str(f"a{i+1}") for i in range(len(args))] + elem_class = Ele(type, type, args=arg_spec) + for tag, nodes in cells.items(): elem = elem_class([], *args, name=tag) # for node in nodes: # print(node, self.get_node(node)) # print("") self.elem(elem, list(map(int,nodes)), tag=tag) + + def apply(self, prototypes=None, **kwds): if prototypes is None: prototypes = self.prototypes From d59b110443a9edce5b640db9f2b16324e96197ee Mon Sep 17 00:00:00 2001 From: claudio perez <50180406+claudioperez@users.noreply.github.com> Date: Tue, 20 Aug 2024 15:43:15 -0700 Subject: [PATCH 6/8] cmp - clean up, organize --- src/opensees/emit/apidoc.py | 37 ++++++++++++------ src/opensees/library/__init__.py | 66 ++------------------------------ src/opensees/library/ast.py | 37 ++++++++++-------- src/opensees/section/__init__.py | 11 +++++- 4 files changed, 61 insertions(+), 90 deletions(-) diff --git a/src/opensees/emit/apidoc.py b/src/opensees/emit/apidoc.py index d590578cdd..d480df9dbd 100644 --- a/src/opensees/emit/apidoc.py +++ b/src/opensees/emit/apidoc.py @@ -4,7 +4,6 @@ import textwrap import importlib -from opensees.emit.writer import ModelWriter from opensees.emit.emitter import Emitter, ScriptBuilder from opensees.library import ast from opensees.library.ast import Arg @@ -13,7 +12,10 @@ 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)+"]" -def write_obj(v, w, qual=None): +def write_obj(v, w, qual=None, ndm=None): + + args = v._args if ndm is None else filter(lambda arg: arg.reqd or ndm in arg.ndm_reqd, v._args) + name = v.__name__ if qual is not None: #qual = qual + "." @@ -21,11 +23,13 @@ def write_obj(v, w, qual=None): else: qual = "" # s = "(" + ", ".join(arg.name for arg in v._args if arg.reqd) + ", **kwds)" - s = "\", " + ", ".join(arg.name for arg in v._args if arg.reqd) + ")" #", **kwds)" + s = "\", " + ", ".join(arg.name for arg in args) + ", **extra)" #s = str(inspect.signature(v)).replace('=None','') + + # if the signature is too long, add in line breaks if len(s) > 45: - s = s.replace(", ", ",
   ") + s = s.replace(", ", ",
    ") w.write(textwrap.dedent(f""" @@ -150,11 +154,18 @@ class ApiDocWriter(ScriptBuilder): def __init__(self): ScriptBuilder.__init__(self, ApiEmitter) - def send(self, obj, idnt=None, qual=None): + def send(self, obj, idnt=None, qual=None, ndm=None): w = self.streams[0] - write_obj(obj, w, qual=qual) - w.endln(); + if ndm is not None: + for dm in ndm: + write_obj(obj, w, qual=qual, ndm=dm) + w.endln(); + + else: + write_obj(obj, w, qual=qual) + w.endln(); + for arg in obj._args: w.write("") @@ -179,12 +190,16 @@ def send(self, obj, idnt=None, qual=None): if __name__ == "__main__": import opensees - _, *module, obj = sys.argv[1].split(".") + _, *module, name = sys.argv[1].split(".") module = ".".join(module) - print(ApiDocWriter().send( - getattr(getattr(opensees, module), obj), qual="Model."+module - )) + obj = getattr(getattr(opensees, module), name) + + ndm = getattr(obj, "_ndms", None) + + print(ApiDocWriter().send(obj, qual="Model."+module, ndm=ndm)) + + diff --git a/src/opensees/library/__init__.py b/src/opensees/library/__init__.py index bbe9494c8e..7515c8b167 100644 --- a/src/opensees/library/__init__.py +++ b/src/opensees/library/__init__.py @@ -25,6 +25,9 @@ def Yld(nature="strength", about=None, **kwds): def Area(**kwds): return Num("A", field="area", about="cross-sectional area", **kwds) +def Izz(**kwds): + return Num("Iz", field="izz", about="area moment of inertia", **kwds) + Rho = lambda: None class ConstitutiveData: _args = [Yng(), Yld(), Rho()] @@ -616,8 +619,6 @@ def create(cls, name, tag, *args): # ]) class element: - Iyc = lambda: Num("iyc", field="iyc", about="Centroidal moment of inertia", alt="section") - Ixc = lambda: Num("ixc", field="ixc", about="", alt="section") ZeroLength = ZeroLength3D = Ele("ZeroLength3D", "zeroLength", @@ -711,7 +712,7 @@ def section(self): def init(self): pass - DisplBeamColumn = dispBeamColumn = Ele("DispBeamColumn", "dispBeamColumn", + DisplBeamColumn = Ele("DispBeamColumn", "dispBeamColumn", #, eleTag, *eleNodes, transfTag, integrationTag, '-cMass', '-mass', mass=0.0) about="Create a dispBeamColumn element.", args=[ @@ -730,65 +731,6 @@ def init(self): refs=["transform"], ) - ElasticBeam2D = Ele("ElasticBeam2D", - "elasticBeamColumn", - args = [ - Tag(), - Grp("nodes", args=[ - Ref("iNode", type=Node, attr="name", about=""), - Ref("jNode", type=Node, attr="name", about=""), - ]), - Area(alt="section"), - Yng( alt="material"), - Iyc(), - Ixc(), - Ref("geom", field="transform", type=Trf, attr="name", default=LinearTransform()), - Num("mass",field="mass_density", flag="-mass", default=0.0, reqd=False, - about="element mass per unit length"), - Flg("-cMass", field="consistent_mass", - about="Flag indicating whether to use consistent mass matrix.") - ], - refs=["transform"], - alts=[ - Ref("material", type="Material"), - Ref("section", type=Sec) - ], - inherit=[_LineElement], - ) - - ElasticBeamColumn3D = Ele("ElasticBeamColumn3D", - "elasticBeamColumn", - args = [ - Tag(), - Grp("nodes", args=[ - Ref("iNode", type=Node, attr="name", about=""), - Ref("jNode", type=Node, attr="name", about=""), - ]), - Area(alt="section"), - Yng( alt="material"), - Num("G", field="shear_modulus", about="", alt="material"), - Num("J", field="torsion_modulus", about="", alt="section"), - #Grp("moi", ctype="struct", args=[ - #Num("iyc", field="iyc", about="Centroidal moment of inertia", alt="section"), - #Num("ixc", field="ixc", about="", alt="section"), - Iyc(), - Ixc(), - #]), - Ref("geom", field="transform", type=Trf, attr="name"), - Num("mass",field="mass_density", flag="-mass", default=0.0, reqd=False, - about="element mass per unit length"), - Flg("-cMass", field="consistent_mass", - about="Flag indicating whether to use consistent mass matrix.") - ], - refs=["transform"], - alts=[ - Ref("material", type="Material"), - Ref("section", type=Sec) - ], - inherit=[_LineElement], - ) - - class constraint: RigidBeamLink = Lnk("RigidBeamLink", diff --git a/src/opensees/library/ast.py b/src/opensees/library/ast.py index eab7f33585..76f99bc51c 100644 --- a/src/opensees/library/ast.py +++ b/src/opensees/library/ast.py @@ -7,28 +7,34 @@ def __repr__(self): return f"${self.name}" class Arg: - __slots__ = ["name", "flag", "value", "field", "default", "type", "reqd", "namespace", "kwds", "about"]#, "group"] + __slots__ = ["name", "flag", "value", "field", "default", "type", + "reqd", "ndm_reqd", "ndm_only", + "namespace", "kwds", "about"]#, "group"] def __init__(self, name = None, #help = None, - flag = None, - reqd = True, - type = None, - field= None, - about= "", - default = None, - # group = None, + flag: str = None, + reqd: bool = True, + ndm_reqd = (), + ndm_only = (), + type = None, + field: str = None, + about = "", + default = None, + # group = None, **kwds ): if name and name[0] == "-": flag = name name = name[1:] - self.name = name - self.flag = flag if flag is not None else "" - self.value = None - self.field = field if field is not None else name - self.type = type - self.reqd = reqd + self.name = name + self.flag = flag if flag is not None else "" + self.value = None + self.field = field if field is not None else name + self.type = type + self.reqd = reqd + self.ndm_reqd = ndm_reqd + self.ndm_only = ndm_only self.default = default self.kwds = kwds self.about = about @@ -65,7 +71,8 @@ def as_tcl_list(self, value=None)->list: return [] else: value = f"${value.name}" - return [self.flag] + [value] + + return [self.flag] + [value] def m_src(self, value=None): value = self._get_value(None,value) diff --git a/src/opensees/section/__init__.py b/src/opensees/section/__init__.py index b28e254a87..216c57ce77 100644 --- a/src/opensees/section/__init__.py +++ b/src/opensees/section/__init__.py @@ -1,13 +1,19 @@ from .section import ( patch, - layer, FiberSection, + ElasticSection +) +layer = patch.layer + +from .shapes import ( ConfinedPolygon, ConfiningPolygon, - PolygonRing + PolygonRing, + sect2gmsh ) SectionGeometry = FiberSection +PlaneShape = FiberSection def _WideFlange(aisc_data, mesh_data, material, tag=None, ndm=None)->SectionGeometry: @@ -70,6 +76,7 @@ def _WideFlange(aisc_data, mesh_data, material, tag=None, ndm=None)->SectionGeom patch.rect(corners=[[ zoff-bi/4,-yoff-tf/2],[ zoff+bi/4, -yoff+tf/2]], material=material, divs=(nfl, nft), rule=int_typ), ]) + def from_shape(type, identifier: str, material=None, mesh=None, units=None, ndm=None, tag=None, **kwds): if identifier == "WF": # TODO From a42de4dfb096bc3dcdb81848932aa78a766618f9 Mon Sep 17 00:00:00 2001 From: claudio perez <50180406+claudioperez@users.noreply.github.com> Date: Tue, 20 Aug 2024 15:43:38 -0700 Subject: [PATCH 7/8] cmp - update OpenSeesRT submodule --- src/libg3 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libg3 b/src/libg3 index 4308e21504..4ded3ac04c 160000 --- a/src/libg3 +++ b/src/libg3 @@ -1 +1 @@ -Subproject commit 4308e215041c66916aee5ce0056cdbfd167b2707 +Subproject commit 4ded3ac04c01bbd84af675f6f9ad5d193a7c8882 From 8c05c45c6f046b365bd324dc9a3d4683951b3d1c Mon Sep 17 00:00:00 2001 From: claudio perez <50180406+claudioperez@users.noreply.github.com> Date: Tue, 20 Aug 2024 15:44:53 -0700 Subject: [PATCH 8/8] cmp --- README.md | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index f823a33307..9b7dee9e05 100644 --- a/README.md +++ b/README.md @@ -30,13 +30,13 @@ Fast nonlinear finite element analysis. `opensees` is a Python package that provides an intuitive API for nonlinear -finite element analysis, implemented in C++. The library features -state-of-the-art finite element formulations and solution algorithms, including -mixed formulations for beams and solids, over 200 material models, and an +finite element analysis, implemented in C++ through the OpenSees framework. +OpenSees features state-of-the-art finite element formulations and solution +algorithms, including mixed formulations for beams and solids, over 200 material models, and an extensive collection of continuation algorithms to solve highly nonlinear problems. -The package supports high quality interactive post processing via the +The `opensees` package supports high quality interactive post processing via the [`sees`](https://pypi.org/project/sees) package. @@ -87,11 +87,21 @@ Additional features include: ### Getting Started +The `opensees` package can be installed into a Python environment +in the standard manner. For example, using `pip`: + +```shell +pip install opensees +``` + +There are several ways to use the `opensees` package: + - To execute Tcl procedures from a Python script, just create an instance - of the `opensees.tcl.Interpreter` class: - ``` + of the `opensees.tcl.Interpreter` class and call its `eval()` method: + ```python interp = opensees.tcl.Interpreter() - interp.eval("model Basic -ndm 2; print -json") + interp.eval("model Basic -ndm 2") + interp.eval("print -json") ``` - To start an interactive interpreter run the shell command: @@ -99,6 +109,10 @@ Additional features include: ```bash python -m opensees ``` + To quit the interpreter, just run `exit`: + ```tcl + opensees > exit + ``` - The `opensees` package exposes a compatibility layer that exactly reproduces the *OpenSeesPy* functions, but does so without mandating a single @@ -115,7 +129,7 @@ Additional features include: `model` function; documentation is under development. -## Compiling +## Development To compile the project see [help/compiling](https://github.com/claudioperez/opensees/blob/master/help/compiling.md)