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)
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
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/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/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
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:
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
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()
-
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)
#
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