"""This module provides a base class for Building Nonlinear
Structural Models (BNSM) to build and export their OpenSees
representations.
"""
# Imports from installed packages
from abc import ABC
import numpy as np
from typing import Optional, List, Dict, Literal, Tuple, Type
import openseespy.opensees as ops
from pathlib import Path
# Imports from bnsm library
from .constants import (RIGID_MAT, BIG_VALUE, GRAV_TS_TAG, GRAV_P_TAG,
NSPA_TS_TAG, NSPA_P_TAG, NSPA_EQUILIBRIUM_TOL,
NSPA_MAX_STIFFNESS_RATIO,
MIN_EIGEN_MODES,
MIN_EIGEN_VALUE)
from .foundation import FoundationBase
from .joint import StairsJointBase, FloorJointBase
from .floor import FloorDiaphragmBase
from .beam import BeamBase, BeamDesign
from .column import ColumnBase, ColumnDesign
from .infill import InfillBase
# Imports from bdim base library
from ...bdim.baselib.building import BuildingBase as BuildingDesign
# Imports from utils library
from ....utils.misc import make_dir, PRECISION
from ....utils import plotter as pl
from ....utils.tags import sub_tag
[docs]
class BuildingBase(ABC):
"""Abstract Base Class for Building Nonlinear Structural Model (BNSM)
implementations. It provides methods to define a BNSM in the OpenSees
domain and to export equivalent Python and Tcl commands.
Attributes
----------
design : ~simdesign.rcmrf.bnsm.baselib.building.BuildingBase
Instance of building design information model.
foundations : List[FoundationBase]
List of foundation instances.
floors : List[FloorDiaphragmBase]
List of floor instances.
floor_joints : List[FloorJointBase]
List of floor joints instances.
stairs_joints : List[StairsJointBase]
List of stairs joints instances.
beams : List[~simdesign.rcmrf.bnsm.baselib.beam.BeamBase]
List of beam instances.
columns : List[~simdesign.rcmrf.bnsm.baselib.column.ColumnBase]
List of column instances.
infills : List[~simdesign.rcmrf.bnsm.baselib.infill.InfillBase]
List of infill instances.
FoundationClass : type[FoundationBase]
Foundation class used to instantiate foundation objects.
FloorClass : type[FloorDiaphragmBase]
Floor diaphragm class used to instantiate floor objects.
FloorJointClass : type[FloorJointBase]
Floor-joint class used to instantiate floor joint objects.
StairsJointClass : type[StairsJointBase]
Stairs-joint class used to instantiate stairs joint objects.
BeamClass : type[~simdesign.rcmrf.bnsm.baselib.beam.BeamBase]
Beam class used by the building factory to instantiate beam objects.
ColumnClass : type[~simdesign.rcmrf.bnsm.baselib.column.ColumnBas]
Column class used by the building factory to instantiate column
objects.
InfillClass : type[~simdesign.rcmrf.bnsm.baselib.infill.InfillBase]
Infill class used to instantiate infill wall objects.
load_factors : Dict[Literal['G', 'Q'], float]
Load factors used to compute gravity loads.
- 'G' : Permanent load factor
- 'Q' : Variable load factor
mass_sources : Dict[Literal['G', 'Q'], float]
Mass sources used to compute seismic masses.
- 'G' : Permanent mass factor
- 'Q' : Variable mass factor
basic_masses : np.ndarray
Point masses for the basic building geometry.
scheme : Literal['FMP', 'EQL', 'MPP', 'TRI', 'UNI']
The loading scheme considered for retriving pushover loads.
- 'FMP : Fundamental-mode proportional loading.
- 'EQL : Equivalent lateral loading.
- 'MPP : Mass proportional loading.
- 'TRI : Triangular or height propotional loading.
- 'UNI : Uniform loading.
max_drift : float
The drift value used to calculate maximum disp. of control node.
dincr: float
First displacement increment considered during nspa.
include_infills: bool
Flag to check whether include infills in the frame model or not.
cyclic_model : bool
If True, the model parameters will be adjusted for cyclic analysis.
cracked_beam : None | bool | float
Cracked-section (effective) flexural properties of the elastic
beam sections: True for My / phi_y of the plastic hinges, or a float
as the ratio of effective to gross stiffness. If None or False,
gross-section properties are used.
cracked_column : None | bool | float
Cracked-section (effective) flexural properties of the elastic
column sections: True for My / phi_y of the plastic hinges, or a
float as the ratio of effective to gross stiffness. If None or
False, gross-section properties are used.
beam_loads : Literal['auto', 'uniform', 'point']
How the beam gravity loads are applied: 'uniform' element loads,
'point' loads at the beam element ends, or 'auto' ('point' if the
beam-column joints are inelastic, 'uniform' otherwise).
"""
design: BuildingDesign
foundations: List[FoundationBase]
floors: List[FloorDiaphragmBase]
floor_joints: List[FloorJointBase]
stairs_joints: List[StairsJointBase]
beams: List[BeamBase]
columns: List[ColumnBase]
infills: List[InfillBase]
FoundationClass: Type[FoundationBase]
FloorDiaphragmClass: Type[FloorDiaphragmBase]
FloorJointClass: Type[FloorJointBase]
StairsJointClass: Type[StairsJointBase]
BeamClass: Type[BeamBase]
ColumnClass: Type[ColumnBase]
InfillClass: Type[InfillBase]
load_factors: Dict[Literal['G', 'Q'], float]
mass_sources: Dict[Literal['G', 'Q'], float]
basic_masses: np.ndarray
scheme: Literal['FMP', 'EQL', 'MPP', 'TRI', 'UNI']
max_drift: float
dincr: float
include_infills: bool
cyclic_model: bool
cracked_beam: Optional[bool | float]
cracked_column: Optional[bool | float]
beam_loads: Literal['auto', 'uniform', 'point']
def __init__(
self, design: BuildingDesign,
load_factors: Dict[Literal['G', 'Q'], float] = {'G': 1.0, 'Q': 0.3},
mass_factors: Dict[Literal['G', 'Q'], float] = {'G': 1.0, 'Q': 0.3},
scheme: Literal['FMP', 'EQL', 'MPP', 'TRI', 'UNI'] = 'EQL',
max_drift: float = 0.05, dincr: float = 0.001,
include_infills: bool = True, cyclic_model: bool = False,
cracked_beam: Optional[bool | float] = None,
cracked_column: Optional[bool | float] = None,
beam_loads: Literal['auto', 'uniform', 'point'] = 'auto',
) -> None:
"""Initialize BNSM object.
Parameters
----------
design : ~simdesign.rcmrf.bnsm.baselib.building.BuildingBase
Instance of building design information model (BDIM)
load_factors : Dict[Literal['G', 'Q'], float], optional
Load factors used to compute gravity loads & seismic masses. By
default {'G': 1.0, 'Q': 0.3}.
- 'G' : Permanent load factor
- 'Q' : Variable load factor
mass_factors : Dict[Literal['G', 'Q'], float], optional
Mass factors used to compute seismic masses. By default
{'G': 1.0, 'Q': 0.3}.
- 'G' : Permanent mass factor
- 'Q' : Variable mass factor
scheme : Literal['FMP', 'EQL', 'MPP', 'TRI', 'UNI'], optional
The loading scheme considered for retriving pushover loads. By
default 'EQL'.
- 'FMP : Fundamental-mode proportional loading.
- 'EQL : Equivalent lateral loading.
- 'MPP : Mass proportional loading.
- 'TRI : Triangular or height propotional loading.
- 'UNI : Uniform loading.
max_drift : float, optional
The drift value used to calculate maximum disp. of control node.
By default 0.05.
dincr: float
First displacement increment considered during nspa.
By default 0.001.
include_infills : bool
Flag to check whether include infills in the frame model or not.
cyclic_model : bool, optional
If True, the model parameters will be adjusted for cyclic analysis.
By default False.
cracked_beam : None | bool | float, optional
Cracked-section (effective) flexural properties of the elastic
beam sections: True for My / phi_y of the plastic hinges, or a
float in (0, 1] as the ratio of effective to gross flexural
stiffness. If None or False, gross-section properties are used.
By default None.
cracked_column : None | bool | float, optional
Cracked-section (effective) flexural properties of the elastic
column sections: True for My / phi_y of the plastic hinges, or a
float in (0, 1] as the ratio of effective to gross flexural
stiffness. If None or False, gross-section properties are used.
By default None.
beam_loads : Literal['auto', 'uniform', 'point'], optional
How the beam gravity loads are applied. By default 'auto'.
- 'uniform' : Uniformly distributed element loads, with the
load over the rigid joint regions applied at the element ends.
- 'point' : Equivalent point loads at the beam element ends
(half of the total beam load at each end), as in O'Reilly and
Sullivan (2019). The beams then carry no gravity moments.
- 'auto' : 'point' if the beam-column joints are modelled as
inelastic, 'uniform' otherwise.
"""
self.design = design
self.load_factors = load_factors
self.mass_sources = mass_factors
self.scheme = scheme
self.max_drift = max_drift
self.dincr = dincr
self.include_infills = include_infills
self.cyclic_model = cyclic_model
self.cracked_beam = self._check_cracked('cracked_beam', cracked_beam)
self.cracked_column = self._check_cracked(
'cracked_column', cracked_column)
self.beam_loads = beam_loads
self.foundations = []
self.floors = []
self.floor_joints = []
self.stairs_joints = []
self.beams = []
self.columns = []
self.infills = []
self._set_basic_masses()
self._initialize_floor_joints()
self._initialize_stairs_joints()
self._initialize_foundations()
self._initialize_beams()
self._initialize_columns()
self._initialize_infills()
@staticmethod
def _check_cracked(
name: str, value: Optional[bool | float]
) -> Optional[bool | float]:
"""Validates a cracked-section input (`cracked_beam` or
`cracked_column`).
Parameters
----------
name : str
Name of the input.
value : None | bool | float
Input value.
Returns
-------
None | bool | float
The value, with numbers converted to float.
Raises
------
ValueError
If the value is not None, a bool or a number in (0, 1].
"""
if value is None or isinstance(value, (bool, np.bool_)):
return None if value is None else bool(value)
if isinstance(value, (int, float, np.integer, np.floating)) and (
0.0 < value <= 1.0):
return float(value)
raise ValueError(
f"Invalid value for '{name}': {value!r}. "
"Use None, a bool or a float in (0, 1].")
@property
def beam_load_type(self) -> Literal['uniform', 'point']:
"""How the beam gravity loads are applied, with 'auto' resolved from
the flexibility model of the floor joints.
Returns
-------
Literal['uniform', 'point']
'point' if ``beam_loads='auto'`` and the joints are inelastic,
'uniform' if ``beam_loads='auto'`` otherwise, else ``beam_loads``.
"""
if self.beam_loads != 'auto':
return self.beam_loads
inelastic = any(joint.flexibility_model == 'inelastic'
for joint in self.floor_joints)
return 'point' if inelastic else 'uniform'
def _initialize_beams(self) -> None:
"""Initialize beam models.
"""
# Initialize the beams
bondslip_fact = self.design.quality.model.bondslip_factor
self.beams = [
self.BeamClass(beam, bondslip_fact, self.load_factors,
self.cyclic_model, self.cracked_beam,
self.beam_load_type)
for beam in self.design.beams
]
# Initialize beam nodes at floor levels. Beams frame into the
# beam-side joint node (floor_node) and columns into center_node, so
# that the joint flexibility element acts as a scissors spring in
# series between them (O'Reilly & Sullivan, 2019). For rigid joints
# floor_node is center_node.
for joint in self.floor_joints:
if joint.design.left_beam:
design = joint.design.left_beam
beam = self._find_beam_by_design(design)
if beam:
beam.ele_node_j = joint.floor_node
beam.jnt_offsets[3] -= joint.bx / 2
if joint.design.right_beam:
design = joint.design.right_beam
beam = self._find_beam_by_design(design)
if beam:
beam.ele_node_i = joint.floor_node
beam.jnt_offsets[0] += joint.bx / 2
if joint.design.rear_beam:
design = joint.design.rear_beam
beam = self._find_beam_by_design(design)
if beam:
beam.ele_node_j = joint.floor_node
beam.jnt_offsets[4] -= joint.by / 2
if joint.design.front_beam:
design = joint.design.front_beam
beam = self._find_beam_by_design(design)
if beam:
beam.ele_node_i = joint.floor_node
beam.jnt_offsets[1] += joint.by / 2
# Initialize beam nodes at mid-storeys (stairs beams)
for joint in self.stairs_joints:
if joint.design.left_beam:
design = joint.design.left_beam
beam = self._find_beam_by_design(design)
if beam:
beam.ele_node_j = joint.center_node
beam.jnt_offsets[3] -= joint.bx / 2
if joint.design.right_beam:
design = joint.design.right_beam
beam = self._find_beam_by_design(design)
if beam:
beam.ele_node_i = joint.center_node
beam.jnt_offsets[0] += joint.bx / 2
def _initialize_columns(self) -> None:
"""Initialize column models.
"""
# Set bondslip factor
bondslip_fact = self.design.quality.model.bondslip_factor
# Check if capacity design is followed for shear
if self.design.OVERSTRENGTH_FACTOR_COLUMN_SHEAR:
capacity_des = True
else:
capacity_des = False
# Initialize columns
self.columns = [
self.ColumnClass(
column, bondslip_fact, capacity_des, self.load_factors,
self.cyclic_model, self.cracked_column
)
for column in self.design.columns
]
# Set rigid offsets at stairs and floor joints
for joint in self.floor_joints + self.stairs_joints:
if joint.design.bottom_column:
design = joint.design.bottom_column
column = self._find_column_by_design(design)
if column:
column.jnt_offsets[5] -= joint.h / 2
column.ele_node_j = joint.center_node
if joint.design.top_column:
design = joint.design.top_column
column = self._find_column_by_design(design)
if column:
column.jnt_offsets[2] += joint.h / 2
column.ele_node_i = joint.center_node
# Set rigid offsets at foundations
for foundation in self.foundations:
if foundation.design.top_column:
design = foundation.design.top_column
column = self._find_column_by_design(design)
if column:
column.ele_node_i = foundation.foundation_node
def _initialize_floor_joints(self) -> None:
"""initialize joint models at floor levels.
"""
# Joint flexibility model
if self.design.beam_type == 1: # Wide beams are used.
# Rigid joints are imposed in the case of wide beams to
# force damage localization.
"""TODO: I am not sure if this should be done here.
Exterior beams are always emergent in fact. We may still see
joint damage. In the future, we may do this check inside, and
determine the modelling option for each joint separately.
Note that the modelling optional is then direction dependent."""
flex_model = 'rigid'
else: # Based on quality and design class for emergent beam cases.
flex_model = self.design.quality.model.joint
floor_i = 0
for points in self.design.geometry.floor_level_points:
floor_i += 1
cnodes = []
joint_masses = []
for point in points:
point_i = self.design.geometry.points.index(point)
mass = self.basic_masses[point_i]
design = self.design._find_joint_by_point(point)
joint = self.FloorJointClass(design, mass, flex_model,
self.load_factors)
cnodes.append(joint.floor_node)
joint_masses.append(mass)
self.floor_joints.append(joint)
# Initiate floor class and append to floors list
self.floors.append(self.FloorDiaphragmClass(
cnodes, floor_i, joint_masses))
def _initialize_stairs_joints(self) -> None:
"""Initialize joint models of stairs (mid-storey level joints).
"""
for point in self.design.geometry.points_at_mid_floor_levels:
point_i = self.design.geometry.points.index(point)
mass = self.basic_masses[point_i]
design = self.design._find_joint_by_point(point)
joint = self.StairsJointClass(design, mass)
self.stairs_joints.append(joint)
def _initialize_foundations(self) -> None:
"""Initialize foundation models.
"""
for point in self.design.geometry.ground_level_points:
point_i = self.design.geometry.points.index(point)
mass = self.basic_masses[point_i]
design = self.design._find_joint_by_point(point)
foundation = self.FoundationClass(design, mass)
self.foundations.append(foundation)
def _set_basic_masses(self) -> None:
"""Sets nodal masses of basic geometry based on mass factors.
"""
elastic_model = self.design.ElasticModelClass(
self.design.beams, self.design.columns, self.design.loads,
self.design.geometry, self.design.beta)
masses_g, masses_q = elastic_model._get_nodal_masses()
self.basic_masses = np.array(masses_g) * self.mass_sources['G'] + \
np.array(masses_q) * self.mass_sources['Q']
def _find_beam_by_design(self, design: BeamDesign) -> Optional[BeamBase]:
"""Finds the beam model by the given design.
Parameters
----------
design : BeamDesign
Beam design instance used for search.
Returns
-------
BeamBase | None
Returns Beam object if design attribute matches with
given design, otherwise, returns None.
"""
filtered_beams = filter(
lambda beam: beam.design == design, self.beams)
# Retrieve the first joint matching the condition
matching_beam = next(filtered_beams, None)
return matching_beam
def _find_column_by_design(
self, design: ColumnDesign
) -> Optional[ColumnBase]:
"""Finds the column model by the given design.
Parameters
----------
design : ColumnDesign
Column design instance used for search.
Returns
-------
ColumnBase | None
Returns Column object if design attribute matches with
given design, otherwise, returns None.
"""
filtered_columns = filter(
lambda column: column.design == design, self.columns)
# Retrieve the first joint matching the condition
matching_column = next(filtered_columns, None)
return matching_column
def _initialize_infills(self) -> None:
"""Initialize infill models.
"""
for design in self.design.infills: # go through each infill design
if isinstance(design.columns[0], list):
c1b = self._find_column_by_design(design.columns[0][0])
c1t = self._find_column_by_design(design.columns[0][1])
s1_ni = c1b.ele_node_i
s2_ni = c1t.ele_node_j
else:
c1 = self._find_column_by_design(design.columns[0])
s1_ni = c1.ele_node_i
s2_ni = c1.ele_node_j
if isinstance(design.columns[1], list):
c2b = self._find_column_by_design(design.columns[1][0])
c2t = self._find_column_by_design(design.columns[1][1])
s1_nj = c2t.ele_node_j
s2_nj = c2b.ele_node_i
else:
c2 = self._find_column_by_design(design.columns[1])
s1_nj = c2.ele_node_j
s2_nj = c2.ele_node_i
strut1_nodes = [s1_ni, s1_nj]
strut2_nodes = [s2_ni, s2_nj]
infill = self.InfillClass(design, strut1_nodes, strut2_nodes)
self.infills.append(infill)
def _add_shared_ops_objects(self) -> None:
"""Add shared OpenSees objects to the domain.
Defines materials, sections, and geometric transformations
that can be referenced by multiple structural components.
"""
# Rigid-like material
ops.uniaxialMaterial('Elastic', RIGID_MAT, BIG_VALUE)
[docs]
def build(self) -> None:
"""Adds the numerical model to the OpenSees domain.
"""
# Destroy any constructed OpenSees object
ops.wipe()
# Create ModelBuilder (with three-dimensions and 6 DOF/node)
ops.model('basic', '-ndm', 3, '-ndf', 6)
# Rigid-like material
self._add_shared_ops_objects()
# Add foundations to ops domain (nodes and support constraints)
for foundation in self.foundations:
foundation.add_to_ops()
# Add stairs joints to ops domain (stair joint nodes)
for joint in self.stairs_joints:
joint.add_to_ops()
# Add floor joints to ops domain (floor cnodes & flexibility elements)
for joint in self.floor_joints:
joint.add_to_ops()
# Add floors to ops domain (floor rnodes & diaphram constraints)
for floor in self.floors:
floor.add_to_ops()
# Add beams to ops domain (beam elements)
for beam in self.beams:
beam.add_to_ops()
# Add columns to ops domain (column elements)
for column in self.columns:
column.add_to_ops()
# Add gravity time-series and load pattern to ops domain
ops.timeSeries('Linear', GRAV_TS_TAG)
ops.pattern('Plain', GRAV_P_TAG, GRAV_TS_TAG)
# Add beam gravity loads to ops domain
for beam in self.beams:
beam.add_grav_loads_to_ops()
# Add column gravity loads to ops domain
for column in self.columns:
column.add_grav_loads_to_ops()
# Perform gravity analysis and save the model state
ops.system('UmfPack')
ops.numberer('RCM')
ops.constraints('Transformation')
ops.test('NormDispIncr', 1e-08, 10, 3)
ops.integrator('LoadControl', 0.1)
ops.algorithm('Newton')
ops.analysis('Static')
ok = ops.analyze(10)
if ok != 0:
raise RuntimeError(
"Gravity analysis did not converge; the model state is "
"invalid for subsequent analyses.")
ops.loadConst('-time', 0.0)
ops.wipeAnalysis()
# Add infills to ops domain (diagonal struts) after gravity
if self.include_infills:
for infill in self.infills:
infill.add_to_ops()
[docs]
def do_modal(
self, num_modes: int = 3, out_dir: Optional[str | Path] = None,
print_screen: bool = False, normalisation: bool = True
) -> Dict[str, List[float]]:
"""Performs modal analysis.
Parameters
----------
num_modes : int, optional
Number of modes for whose properties are calculated. By default 3.
out_dir : bool, optional
Output directory to save the modal properties and eigen vectors
for the floor retained nodes. If None, these are not saved.
By default None.
print_screen : bool, optional
Flag to print modal properties on the screen. By default False.
normalisation : bool, optional
Flag to compute the modal properties by using a
displacement-normalized version of the eigenvectors.
By default True.
Returns
-------
Dict[str, List[float]]
Dictionary containing modal properties. The per-mode entries
(``eigen*`` and ``parti*``) are given for the first `num_modes`
modes.
Notes
-----
At least ``MIN_EIGEN_MODES`` modes are solved for, because the
default ARPACK solver may not converge for a single eigenvalue and
then returns unset values (e.g. periods of 1e157 s). Eigenvalues
below ``MIN_EIGEN_VALUE`` are taken as unset, and the next solver is
tried. The modal properties file written to `out_dir` contains all
solved modes.
TODO
----
For now print screen does not work.
"""
# Build the model
self.build()
# Perform eigen value analysis, with at least MIN_EIGEN_MODES modes
num_solve = max(num_modes, MIN_EIGEN_MODES)
list_solvers = ['-genBandArpack', '-fullGenLapack', '-symmBandLapack']
ok = False
for s in list_solvers:
try:
eigen_values = ops.eigen(s, num_solve)
# All the eigenvalues must be set (not tiny) and finite
ok = all(MIN_EIGEN_VALUE < value < 1e+300
for value in eigen_values)
if ok:
break
except BaseException as e:
print(f"Using {s[1:]} as solver... received an error: {e}")
if not ok:
raise RuntimeError(
"Could not complete the eigenvalue analysis, "
+ "something is wrong...\n"
+ "Try to reduce number of modes to determine..."
)
# Save eigen vectors for retained floor nodes
if out_dir:
# Create the output directory
out_dir = Path(out_dir)
if not Path.exists(out_dir):
make_dir(out_dir)
# Modal properties directory
modal_path = (out_dir / 'ModalProperties.txt').as_posix()
# Save eigen vectors for retained floor nodes
nodes = [floor.rnode.tag for floor in self.floors]
for i in range(1, num_modes + 1):
modal_disps = []
for node in nodes:
disps = ', '.join([
f'{disp}' for disp in ops.nodeEigenvector(node, i)
])
modal_disps.append(f'{node}, {disps}')
modal_disps = '\n'.join(modal_disps)
eigen_path = (out_dir / f'EigenVectors_Mode{i}.txt').as_posix()
with open(eigen_path, 'w') as file:
file.write(modal_disps)
# Set modal analysis arguments
args = []
if print_screen:
args.append('-print')
if out_dir:
args.extend(['-file', modal_path])
if normalisation:
args.append('-unorm')
args.append('-return')
# Perform modal analysis
modal_properties = ops.modalProperties(*args)
# Keep the requested modes in the per-mode properties
for key, values in modal_properties.items():
if key.startswith(('eigen', 'parti')):
modal_properties[key] = values[:num_modes]
return modal_properties
[docs]
def do_nspa(
self, ctrl_dof: Literal[1, 2],
out_dir: Optional[str | Path] = None,
use_recorder: bool = False
) -> Tuple[List[float], List[float], int]:
"""Performs nonlinear static pushover analysis (NSPA).
Parameters
----------
ctrl_dof : Literal[1, 2]
Control degrees of freedom for loading.
1: X-direction.
2: Y-direction.
out_dir : str | Path, optional
Directory to save outputs. If None, no output is saved.
By default None.
use_recorder : bool, optional
If True, OpenSees recorders are used. By default False.
Return
------
ctrl_disp : List[float]
Displacement values of control node.
base_shear : List[float]
Base shear value obtained as sum of the reaction forces.
ok : int
0 if the analysis ends at the target displacement
(``max_drift``) or after a loss of strength, i.e., a converged
step with a base shear below 60% of the maximum (recorded).
Negative if it stops earlier (see Notes).
Notes
-----
A step that fails is repeated with alternative algorithms, smaller
displacement increments and, as a last resort, relaxed tolerances.
The analysis stops before the target displacement, with ``ok``
negative, if no step converges, if a step is out of equilibrium
(base shear differing from the applied load by more than
``NSPA_EQUILIBRIUM_TOL``), if the tangent stiffness of a step
exceeds ``NSPA_MAX_STIFFNESS_RATIO`` times the initial stiffness,
or if the base shear reaches 1.5 times the seismic weight. The
exported ``nspa.py`` and ``nspa.tcl`` do not include the
equilibrium and stiffness checks.
"""
# Analysis direction
direction = 'x' if ctrl_dof == 1 else 'y'
# Output directory is provided, save stuff
if out_dir:
out_dir = Path(out_dir)
# Create the output directory if needed
if not Path.exists(out_dir):
make_dir(out_dir)
# Save the log without any echo
ops.logFile(str(out_dir / f'log_{direction}.out'), '-noEcho')
# Get NSPA loading parameters
nodes, loads, ctrl_node = self._get_nspa_loading_parameters(ctrl_dof)
# Add NSPA time-series and load pattern to ops domain
ops.timeSeries('Linear', NSPA_TS_TAG)
ops.pattern('Plain', NSPA_P_TAG, NSPA_TS_TAG)
# Add lateral nspa loads to ops domain
for node, load_values in zip(nodes, loads):
ops.load(node, *load_values)
# Total lateral load of the pattern (load factor of 1)
total_load = abs(sum(values[ctrl_dof - 1] for values in loads))
# Set foundation and floor nodes
supports = [found.foundation_node.tag for found in self.foundations]
floors = [floor.rnode.tag for floor in self.floors]
# Get the base level and storey heights
base_level = min([ops.nodeCoord(node, 3) for node in supports])
storey_heights = [round(ops.nodeCoord(node, 3) - base_level,
PRECISION) for node in floors]
# Output directory is provided, save stuff
if out_dir:
# Save the storey heights
np.savetxt(out_dir / 'storey_heights.out',
storey_heights, fmt="%.3f")
# Set the recorder for storey displacements and support reactions
if use_recorder:
ops.recorder(
'Node', '-file',
str(out_dir / f'storey_displacements_{direction}.out'),
'-node', *floors, '-dof', ctrl_dof, 'disp')
ops.recorder(
'Node', '-file',
str(out_dir / f'support_reactions_{direction}.out'),
'-node', *supports, '-dof', ctrl_dof, 'reaction')
# Set some analysis parameters
factors = [1.0, 0.5, 0.25, 0.1, 0.01, 0.001]
dincr_values = [factor * self.dincr for factor in factors]
dincr_values.sort(reverse=True) # should be in descending order
dincr = max(dincr_values)
max_disp = self.max_drift * (ops.nodeCoord(ctrl_node, 3) - base_level)
tol_init = 1.0e-6
iter_init = 20
shear_lim = sum(sum(floor.masses) for floor in self.floors) * 10 * 1.5
ops.wipeAnalysis()
ops.system('UmfPack')
ops.numberer('RCM')
ops.constraints('Transformation')
ops.test('EnergyIncr', tol_init, iter_init)
ops.integrator('DisplacementControl', ctrl_node, ctrl_dof, dincr)
ops.algorithm('Newton', '-initialThenCurrent')
ops.analysis('Static')
# Start performing the analysis
ctrl_disp = [0]
base_shear = [0]
storey_displacements = [[0] * len(floors)]
support_reactions = [[0] * len(supports)]
ok = 0
cont = True
relaxed = False
init_stiffness = np.inf # initial stiffness, set after the 1st step
while ok == 0 and cont:
if relaxed: # keep the last algorithm, reset the tolerance
relaxed = False
ops.test('NormDispIncr', tol_init, 100)
# Index for the current displacement increment
dincr_idx = dincr_values.index(dincr)
# Gradually increase the displacement increment size
if dincr_idx != 0:
dincr_idx -= 1
dincr = dincr_values[dincr_idx]
# Set current step size and try performing analysis with it
ops.integrator(
'DisplacementControl', ctrl_node, ctrl_dof, dincr
)
ok = ops.analyze(1)
# Try same dincr with alternative algorithms first
if ok != 0:
ok = self._set_nspa_algorithm(
tol_init, ctrl_node, ctrl_dof, dincr
)
# If still failing, keep decreasing dincr using dincr_values
while ok != 0 and dincr_idx < len(dincr_values) - 1:
dincr_idx += 1
dincr = dincr_values[dincr_idx]
ok = self._set_nspa_algorithm(
tol_init, ctrl_node, ctrl_dof, dincr
)
# If still failing, relax tolerance
if ok != 0: # Increase tolerance by factor of 10
relaxed = True
ok = self._set_nspa_algorithm(
10 * tol_init, ctrl_node, ctrl_dof, dincr, iter=100,
)
if ok != 0: # increase tolerance by factor of 100
ok = self._set_nspa_algorithm(
100 * tol_init, ctrl_node, ctrl_dof, dincr, iter=200
)
# Get the base shear force
ops.reactions()
current_ctrl_disp = ops.nodeDisp(ctrl_node, ctrl_dof)
current_reactions = [
ops.nodeReaction(node, ctrl_dof) for node in supports
]
current_base_shear = abs(sum(current_reactions))
# Stop if a step is accepted out of equilibrium, i.e., the base
# shear differs from the applied load. This can happen with the
# relaxed tolerances, and the analysis would continue from an
# invalid committed state (e.g., unrealistic base shear spikes).
applied_load = abs(ops.getLoadFactor(NSPA_P_TAG)) * total_load
if ok == 0 and abs(current_base_shear - applied_load) > (
NSPA_EQUILIBRIUM_TOL * max(applied_load, 1.0)):
ok = -1
# Check the base shear limits
within_limits = (current_base_shear < shear_lim and
current_base_shear >= 0.6 * max(base_shear))
# Set continue flag
cont = within_limits and current_ctrl_disp < max_disp
# Stop if the tangent stiffness of the step clearly exceeds the
# initial stiffness, which a valid pushover curve never does. Such
# a jump indicates a degenerate element state (unrealistic spikes).
if ok == 0 and len(base_shear) > 1:
dU = abs(current_ctrl_disp - ctrl_disp[-1])
dV = current_base_shear - base_shear[-1]
if dV > NSPA_MAX_STIFFNESS_RATIO * init_stiffness * dU:
ok = -1
# Initial stiffness, from the first step
if ok == 0 and len(base_shear) == 1 and current_ctrl_disp != 0:
init_stiffness = current_base_shear / abs(current_ctrl_disp)
# Append the results, including the step reaching max_disp and
# a converged step where the base shear drops below 60% of the
# maximum (sudden loss of strength, which ends the analysis)
if ok == 0 and current_base_shear < shear_lim:
ctrl_disp.append(current_ctrl_disp)
base_shear.append(current_base_shear)
storey_displacements.append(
[ops.nodeDisp(node, ctrl_dof) for node in floors]
)
support_reactions.append(current_reactions)
elif (
max(base_shear) == 0.0 # Analysis did not even start
or base_shear[-1] / max(base_shear) > 0.8 # Not good enough
):
ok = -1
if out_dir: # Output directory is provided, save stuff
np.savetxt(out_dir / f'ctrl_disp_{direction}.out',
ctrl_disp, delimiter=' ', fmt="%.8f")
np.savetxt(out_dir / f'base_shear_{direction}.out',
base_shear, delimiter=' ', fmt="%.8f")
if not use_recorder: # Recorders are not used
np.savetxt(out_dir / f'storey_displacements_{direction}.out',
storey_displacements, delimiter=' ', fmt="%.8f")
np.savetxt(out_dir / f'support_reactions_{direction}.out',
support_reactions, delimiter=' ', fmt="%.8f")
# Wipe the numerical model
ops.wipe()
# Return base shear and control node displacement history
return ctrl_disp, base_shear, ok
def _set_nspa_algorithm(
self, tol: float, ctrl_node: int, ctrl_dof: int, dincr: float,
iter: int = None
) -> None:
"""Sets the solution algorithm for NSPA in ops domain.
Parameters
----------
tol : float
The tolerance criteria used to check for convergence.
ctrl_node : int
Tag of control node.
ctrl_dof : int
Tag of control degrees of freedom, i.e., 1 or 2.
dincr : float
The displacement increment considered during analysis.
iter : int, optional
The max number of iterations to check before returning failure.
By default 100.
Return
------
int
Result of the new analysis step in OpenSees.
"""
# Per-algorithm iteration budgets
iter_krylov = iter if iter else 50
iter_ls = iter if iter else 50
iter_bfgs = iter if iter else 100
# 1. KrylovNewton
ops.test('NormDispIncr', tol, iter_krylov)
ops.integrator('DisplacementControl', ctrl_node, ctrl_dof, dincr)
ops.algorithm('KrylovNewton')
ok = ops.analyze(1)
# 2. NewtonLineSearch
if ok != 0 and iter_ls >= 50:
ops.test('NormDispIncr', tol, iter_ls)
ops.algorithm('NewtonLineSearch', '-InitialInterpolated', 0.8)
ok = ops.analyze(1)
# 3. BFGS
if ok != 0 and iter_bfgs >= 100:
ops.test('NormDispIncr', tol, iter_bfgs)
ops.algorithm('BFGS')
ok = ops.analyze(1)
# Return the analysis result
return ok
def _get_nspa_loading_parameters(
self, ctrl_dof: Literal[1, 2]
) -> Tuple[List[int], List[float], int]:
"""Gets the parameters required for nspa loading.
Parameters
----------
ctrl_dof : Literal[1, 2]
Control degrees of freedom for loading.
1: X-direction.
2: Y-direction.
Returns
-------
nodes : List[int]
List containing tags of nodes which will be pushed.
loads : List[float]
List of loads which corresponds to each node deemed to be pushed.
ctrl_node : int
Tag of control node (determined based on max height).
"""
# Find all required metrics for load computations.
nodes = []
masses = []
heights = []
modal_disps = []
if self.scheme == 'FMP': # Do Modal Analysis
num_modes = round(len(self.basic_masses) / 2)
modal_properties = self.do_modal(num_modes)
if ctrl_dof == 1:
mass_part_rat = np.array(modal_properties['partiMassRatiosMX'])
elif ctrl_dof == 2:
mass_part_rat = np.array(modal_properties['partiMassRatiosMY'])
inds = mass_part_rat.argsort()[::-1]
mode1 = int(inds[0] + 1)
else: # Just build the model
self.build()
# Base level, floor heights are measured from it
base_level = min(ops.nodeCoord(foundation.foundation_node.tag, 3)
for foundation in self.foundations)
for floor in self.floors:
node = floor.rnode.tag
nodes.append(node)
masses.append(sum(floor.masses))
heights.append(ops.nodeCoord(node, 3) - base_level)
if self.scheme == 'FMP':
modal_disps.append(ops.nodeEigenvector(node, mode1, ctrl_dof))
# Convert these lists to arrays for performing array operations
masses = np.array(masses)
heights = np.array(heights)
modal_disps = np.array(modal_disps)
# Loads for fundamental-mode proportional loading scheme
if self.scheme == 'FMP':
loads = (masses * modal_disps) / np.sum(masses * modal_disps)
# Loads for equivalent lateral loading scheme
elif self.scheme == 'EQL':
loads = (masses * heights) / np.sum(masses * heights)
# Loads for mass proportional loading scheme
elif self.scheme == 'MPP':
loads = masses / np.sum(masses)
# Loads for triangular or height propotional loading scheme
elif self.scheme == 'TRI':
loads = heights / np.sum(heights)
# Loads for uniform loading scheme
elif self.scheme == 'UNI':
loads = np.ones(len(nodes)) / len(nodes)
# Make loading positive
if np.sum(loads) < 0:
loads = -1.0 * loads
# Get control node
inds = heights.argsort()[::-1]
ctrl_node = nodes[inds[0]]
loads = loads.tolist()
if ctrl_dof == 1:
loads = [[load, 0, 0, 0, 0, 0] for load in loads]
elif ctrl_dof == 2:
loads = [[0, load, 0, 0, 0, 0] for load in loads]
return nodes, loads, ctrl_node
[docs]
def to_py(self, directory: Path | str) -> None:
"""Exports the numerical model files in Python format.
Parameters
----------
directory : Path | str
Directory in which the files will be written.
"""
# Create output directory
directory = Path(directory)
if not Path.exists(directory):
make_dir(directory)
# Get each function as list of strings
build_list = self._get_build_py()
foundations_list = self._get_foundations_py()
joints_list = self._get_joints_py()
floors_list = self._get_floors_py()
beams_list = self._get_beams_py()
columns_list = self._get_columns_py()
infills_list = self._get_infills_py()
gravity_list = self._get_gravity_py()
modal_list = self._get_modal_py()
nspa_list = self._get_nspa_py()
run_list = self._get_run_py()
# Write __init__ file
with open(directory / '__init__.py', 'w') as file:
file.write('\n')
# Write model source file
with open(directory / 'model.py', 'w') as file:
file.write('\n'.join(build_list))
# Write foundations source file
with open(directory / 'foundations.py', 'w') as file:
file.write('\n'.join(foundations_list))
# Write joints source file
with open(directory / 'joints.py', 'w') as file:
file.write('\n'.join(joints_list))
# Write floors source file
with open(directory / 'floors.py', 'w') as file:
file.write('\n'.join(floors_list))
# Write beams source file
with open(directory / 'beams.py', 'w') as file:
file.write('\n'.join(beams_list))
# Write columns source file
with open(directory / 'columns.py', 'w') as file:
file.write('\n'.join(columns_list))
# Write infills source file
with open(directory / 'infills.py', 'w') as file:
file.write('\n'.join(infills_list))
# Write gravity analysis source file
with open(directory / 'gravity.py', 'w') as file:
file.write('\n'.join(gravity_list))
# Write modal analysis source file
with open(directory / 'modal.py', 'w') as file:
file.write('\n'.join(modal_list))
# Write nspa source file
with open(directory / 'nspa.py', 'w') as file:
file.write('\n'.join(nspa_list))
# Write example run file
with open(directory / 'run.py', 'w') as file:
file.write('\n'.join(run_list))
[docs]
def to_tcl(self, directory: Path | str) -> None:
"""Exports the numerical model files in Tcl format.
Parameters
----------
directory : Path | str
Directory in which the files will be written.
"""
# Create output directory
directory = Path(directory)
if not Path.exists(directory):
make_dir(directory)
# Get each function as list of strings
build_list = self._get_build_tcl()
foundations_list = self._get_foundations_tcl()
joints_list = self._get_joints_tcl()
floors_list = self._get_floors_tcl()
beams_list = self._get_beams_tcl()
columns_list = self._get_columns_tcl()
infills_list = self._get_infills_tcl()
gravity_list = self._get_gravity_tcl()
modal_list = self._get_modal_tcl()
nspa_list = self._get_nspa_tcl()
run_list = self._get_run_tcl()
# Write model source file
with open(directory / 'model.tcl', 'w') as file:
file.write('\n'.join(build_list))
# Write foundations source file
with open(directory / 'foundations.tcl', 'w') as file:
file.write('\n'.join(foundations_list))
# Write joints source file
with open(directory / 'joints.tcl', 'w') as file:
file.write('\n'.join(joints_list))
# Write floors source file
with open(directory / 'floors.tcl', 'w') as file:
file.write('\n'.join(floors_list))
# Write beams source file
with open(directory / 'beams.tcl', 'w') as file:
file.write('\n'.join(beams_list))
# Write columns source file
with open(directory / 'columns.tcl', 'w') as file:
file.write('\n'.join(columns_list))
# Write infills source file
with open(directory / 'infills.tcl', 'w') as file:
file.write('\n'.join(infills_list))
# Write gravity analysis source file
with open(directory / 'gravity.tcl', 'w') as file:
file.write('\n'.join(gravity_list))
# Write modal analysis source file
with open(directory / 'modal.tcl', 'w') as file:
file.write('\n'.join(modal_list))
# Write nspa source file
with open(directory / 'nspa.tcl', 'w') as file:
file.write('\n'.join(nspa_list))
# Write example run file
with open(directory / 'run.tcl', 'w') as file:
file.write('\n'.join(run_list))
def _get_shared_opspy(self) -> List[str]:
"""Generate Python lines for shared OpenSees domain objects.
Return
------
List[str]
List of lines that define reusable materials, sections, and
geometric transformations that can be referenced by multiple
structural components.
"""
return [
"# Rigid-like material",
f"ops.uniaxialMaterial('Elastic', {RIGID_MAT}, {BIG_VALUE})"
]
def _get_build_py(self) -> List[str]:
"""Generates a list of strings representing the Python function that
constructs the components of the numerical model in OpenSees.
Return
------
List[str]
List of lines which will be written into model.py file.
"""
# Get lines of adding the shared objects
shared = self._get_shared_opspy()
# Add lines for setting build method
content = [
'"""Adds the numerical model to the OpenSees domain.',
'"""',
"ops.wipe()",
"ops.model('basic', '-ndm', 3, '-ndf', 6)",
"",
*shared,
"",
"# Define components of foundations",
"add_foundations()",
"# Define components of joints",
"add_joints()",
"# Define components of floors",
"add_floors()",
"# Define components of beams",
"add_beams()",
"# Define components of columns",
"add_columns()",
"# Perform static analysis under gravity loads",
"do_gravity()",
"# Define components of infills (after gravity, unstressed)",
"if infills:",
" add_infills()",
""
]
# Add white spaces for method content
content = [' ' + item if item
else item
for item in content]
# Add method definition
method = [f"def build_model(infills={self.include_infills}) -> None:"]
# Add imports
imports = ["import openseespy.opensees as ops",
"",
"from .foundations import add_foundations",
"from .floors import add_floors",
"from .joints import add_joints",
"from .beams import add_beams",
"from .columns import add_columns",
"from .infills import add_infills",
"from .gravity import do_gravity",
"", ""]
content = imports + method + content
return content
def _get_foundations_py(self) -> List[str]:
"""Generates a list of strings representing the Python function that
constructs the components of the foundation in OpenSees.
Return
------
List[str]
List of lines which will be written into foundations.py file.
"""
# Initialize the content list - docstring
content = [
'"""Add foundation components to ops domain '
'(nodes and constraints).',
'"""'
]
# Add lines for constructing the foundation objects
for foundation in self.foundations:
content.extend(foundation.to_py())
content.append("")
# Add white spaces for method content
content = [' ' + item if item
else item
for item in content]
# Add method definition
method = ["def add_foundations() -> None:"]
# Add imports
imports = ["import openseespy.opensees as ops",
"", ""]
content = imports + method + content
return content
def _get_floors_py(self) -> List[str]:
"""Generates a list of strings representing the Python function that
constructs the components of all floors in OpenSees.
Return
------
List[str]
List of lines which will be written into floors.py file.
"""
# Initialize the content list - docstring
content = [
'"""Add floors to ops domain (retained nodes & diaphrams).',
'"""'
]
# Add lines for constructing the floor objects
for floor in self.floors:
content.extend(floor.to_py())
content.append("")
# Add white spaces for method content
content = [' ' + item if item
else item
for item in content]
# Add method definition
method = ["def add_floors() -> None:"]
# Add imports
imports = ["import openseespy.opensees as ops",
"", ""]
content = imports + method + content
return content
def _get_joints_py(self) -> List[str]:
"""Generates a list of strings representing the Python function that
constructs the components of all joints in OpenSees.
Return
------
List[str]
List of lines which will be written into joints.py file.
"""
# Initialize the content list - docstring
content = [
'"""Add components of joints to ops domain.',
'"""'
]
# Add lines for constructing the stairs joint objects
content.append('# -------------------------------------------------')
content.append('# Add stairs joints to ops domain')
content.append('# -------------------------------------------------')
for joint in self.stairs_joints:
content.extend(joint.to_py())
content.append("")
# Add lines for constructing the floor joint objects
content.append('# -------------------------------------------------')
content.append('# Add floor joints to ops domain')
content.append('# -------------------------------------------------')
for joint in self.floor_joints:
content.extend(joint.to_py())
content.append("")
# Add white spaces for method content
content = [' ' + item if item
else item
for item in content]
# Add method definition
method = ["def add_joints() -> None:"]
# Add imports
imports = ["import openseespy.opensees as ops",
"", ""]
content = imports + method + content
return content
def _get_beams_py(self) -> List[str]:
"""Generates a list of strings representing the Python function that
constructs the components of all beams in OpenSees.
Return
------
List[str]
List of lines which will be written into beams.py file.
"""
# Initialize the content list - docstring
content = [
'"""Add components of all beams to ops domain',
'"""'
]
# Add lines for constructing the beam objects
for beam in self.beams:
content.extend(beam.to_py())
content.append("")
# Add white spaces for method content
content = [' ' + item if item
else item
for item in content]
# Add method definition
method = ["def add_beams() -> None:"]
# Add imports
imports = ["import openseespy.opensees as ops",
"", ""]
content = imports + method + content
return content
def _get_columns_py(self) -> List[str]:
"""Generates a list of strings representing the Python function that
constructs the components of all columns in OpenSees.
Return
------
List[str]
List of lines which will be written into columns.py file.
"""
# Initialize the content list - docstring
content = [
'"""Add components of all columns to ops domain',
'"""'
]
# Add lines for constructing the beam objects
for column in self.columns:
content.extend(column.to_py())
content.append("")
# Add white spaces for method content
content = [' ' + item if item
else item
for item in content]
# Add method definition
method = ["def add_columns() -> None:"]
# Add imports
imports = ["import openseespy.opensees as ops",
"", ""]
content = imports + method + content
return content
def _get_infills_py(self) -> List[str]:
"""Generates a list of strings representing the Python function that
constructs the components of all the masonry infill walls in OpenSees.
Return
------
List[str]
List of lines which will be written into infills.py file.
"""
# Initialize the content list - docstring
content = [
'"""Add components of all infills to ops domain',
'"""'
]
# Add lines for constructing the beam objects
for infill in self.infills:
content.extend(infill.to_py())
content.append("")
# Add white spaces for method content
content = [' ' + item if item
else item
for item in content]
# Add method definition
method = ["def add_infills() -> None:"]
# Add imports
imports = ["import openseespy.opensees as ops",
"", ""]
content = imports + method + content
return content
def _get_gravity_py(self) -> List[str]:
"""Generates a list of strings representing the Python function that
constructs the components of the gravity analysis in OpenSees.
Return
------
List[str]
List of lines which will be written into gravity.py file.
"""
# Initialize the content list - docstring, load pattern, time-series
content = [
'"""Perform linear static analysis under gravity loads.',
'"""',
"# Add gravity time-series and load pattern to ops domain",
f"ops.timeSeries('Linear', {GRAV_TS_TAG})",
f"ops.pattern('Plain', {GRAV_P_TAG}, {GRAV_TS_TAG})",
""
]
# Add lines for constructing the beam gravity load objects
content.append("# Add beam gravity loads to ops domain")
for beam in self.beams:
content.extend(beam.to_py_grav_loads())
# Add lines for constructing the column gravity load objects
content.append("")
content.append("# Add column gravity loads to ops domain")
for column in self.columns:
content.extend(column.to_py_grav_loads())
# Add lines for performing gravity analysis
content.append("")
content.append("# Perform gravity analysis and save the model state")
content.append("ops.system('UmfPack')")
content.append("ops.numberer('RCM')")
content.append("ops.constraints('Transformation')")
content.append("ops.test('NormDispIncr', 1e-08, 10, 3)")
content.append("ops.integrator('LoadControl', 0.1)")
content.append("ops.algorithm('Newton')")
content.append("ops.analysis('Static')")
content.append("ok = ops.analyze(10)")
content.append("if ok != 0:")
content.append(" raise RuntimeError('Gravity analysis did not "
"converge; the model state is invalid for subsequent "
"analyses.')")
content.append("ops.loadConst('-time', 0.0)")
content.append("ops.wipeAnalysis()")
content.append("")
# Add white spaces for method content
content = [' ' + item if item
else item
for item in content]
# Add method definition
method = ["def do_gravity() -> None:"]
# Add imports
imports = ["import openseespy.opensees as ops",
"", ""]
content = imports + method + content
return content
def _get_modal_py(self) -> List[str]:
"""Generates a list of strings representing the Python function
for modal analysis in OpenSees.
Return
------
List[str]
List of lines which will be written into modal.py file.
"""
# Retrained floor nodes
nodes = ', '.join([f"{floor.rnode.tag}" for floor in self.floors])
# Initialize the content list - docstring
content = ['"""Perform modal analysis for built OpenSees model.',
"",
"Parameters",
"----------",
"num_modes : int, optional",
" Number of modes considered for modal analysis.",
" By default 3.",
"",
"Return",
"------",
"dict",
" Dictionary containing modal properties.",
'"""',]
# Add lines for setting the output directory
content.append("# Set output directory")
content.append(
"output_directory = Path(__file__).parent / 'Modal-Results'")
content.append("if not Path.exists(output_directory):")
content.append(" Path.mkdir(output_directory)")
# Add lines for building the numerical model
content.append("# Build the numerical model")
content.append("build_model()")
# Add lines for performing eigenvalue analysis
content.append("")
content.append("# Perform eigen value analysis, with at least "
f"{MIN_EIGEN_MODES} modes (ARPACK is unreliable for a "
"single eigenvalue)")
content.append(f"eigen_values = ops.eigen(max(num_modes, "
f"{MIN_EIGEN_MODES}))")
content.append(f"if not all({MIN_EIGEN_VALUE} < value < 1e+300 "
"for value in eigen_values):")
content.append(" raise RuntimeError(f'Invalid eigenvalues: "
"{eigen_values}')")
# Add lines for looping through eigenvectors and saving them
content.append("# Save eigen vectors for retained floor nodes")
content.append(f"nodes = [{nodes}]")
content.append("for i in range(1, num_modes+1):")
content.append(" modal_disps = []")
content.append(" for node in nodes:")
content.append(
" disps = "
"', '.join([f'{disp}' for disp in ops.nodeEigenvector(node, i)])"
)
content.append(" modal_disps.append(f'{node}, {disps}')")
content.append(" modal_disps = '\\n'.join(modal_disps)")
content.append(" report_file_path = (output_directory / "
"f'EigenVectors_Mode{i}.txt').as_posix()")
content.append(" with open(report_file_path, 'w') as file:")
content.append(" file.write(modal_disps)")
# Add lines for saving modal properties
content.append("")
content.append("# Perform modal analysis and save results")
content.append("report_file_path = (output_directory / "
"'ModalProperties.txt').as_posix()")
content.append("results = ops.modalProperties('-print', '-return', "
"'-file', report_file_path, '-unorm')")
content.append("# Keep the requested modes in the per-mode properties")
content.append("for key, values in results.items():")
content.append(" if key.startswith(('eigen', 'parti')):")
content.append(" results[key] = values[:num_modes]")
content.append("")
# Add return statement
content.append("return results")
content.append("")
# Add white spaces for method content
content = [' ' + item if item
else item
for item in content]
# Add method definition
method = ["def do_modal(num_modes: int = 3) -> dict:"]
# Add imports
imports = ["import openseespy.opensees as ops",
"from pathlib import Path",
"",
"from .model import build_model",
"", ""]
content = imports + method + content
return content
def _get_nspa_dof_py(self, ctrl_dof: Literal[1, 2]) -> List[str]:
"""Generates a list of strings representing the Python function
for nonlinear static pushover (NSPA) analysis in OpenSees.
Parameters
----------
ctrl_dof : Literal[1, 2]
Control degrees of freedom for loading.
1: X-direction.
2: Y-direction.
Return
------
List[str]
List of lines representing do_nspa_{direction} method.
"""
# Get NSPA loading parameters
nodes, loads, ctrl_node = \
self._get_nspa_loading_parameters(ctrl_dof)
# Set support nodes
supports = ', '.join([
f"{found.foundation_node.tag}" for found in self.foundations
])
# Set floor nodes
floors = ', '.join([f"{floor.rnode.tag}" for floor in self.floors])
# Directions per dof
if ctrl_dof == 1:
direction = 'x'
elif ctrl_dof == 2:
direction = 'y'
# Initialize the content list - docstring
content = [
'"""Performs nonlinear static pushover analysis (NSPA) '
f'in {direction} direction.',
"",
"Parameters",
"----------",
"max_drift : float, optional.",
" Maximum considered drift value for the control node.",
f" By default {self.max_drift}",
"dincr : float, optional.",
" First displacement increment considered during the analysis.",
f" By default {self.dincr}.",
"",
"Return",
"------",
"ctrl_disp : List[float]",
" Displacement values of control node.",
"base_shear : List[float]",
" Base shear value obtained as sum of the reaction forces.",
'"""'
]
# Add lines for setting the output directory
content.append("# Set output directory")
content.append(
"output_directory = Path(__file__).parent / 'NSPA-Results'")
content.append("if not Path.exists(output_directory):")
content.append(" Path.mkdir(output_directory)")
content.append("reaction_file_path = (output_directory / "
f"'support_reactions_{direction}.out').as_posix()")
content.append("disp_file_path = (output_directory / "
f"'storey_displacements_{direction}.out').as_posix()")
content.append("storey_heights_file_path = (output_directory / "
"'storey_heights.out').as_posix()")
# Add lines for building the numerical model
content.append("")
content.append("# Build the numerical model")
content.append("build_model()")
# Add lines for setting the time-series and load pattern
content.append("")
content.append("# Add NSPA time-series and load pattern to ops domain")
content.append(f"ops.timeSeries('Linear', {NSPA_TS_TAG})")
content.append(f"ops.pattern('Plain', {NSPA_P_TAG}, {NSPA_TS_TAG})")
# Add lines for setting the nspa loads
content.append("# Add lateral nspa loads to ops domain")
for node, load_values in zip(nodes, loads):
values_str = ', '.join([f"{val}" for val in load_values])
content.append(f"ops.load({node}, {values_str})")
# Add lines for setting recorders
content.append("")
content.append("# Set the recorders")
content.append(f"ctrl_node = {ctrl_node} # Control node")
content.append(f"ctrl_dof = {ctrl_dof} # Control dof")
content.append(f"supports = [{supports}] # Foundation nodes")
content.append(f"floors = [{floors}] # Retained floor nodes")
content.append("ops.recorder('Node', '-file', disp_file_path, '-node',"
" *floors, '-dof', ctrl_dof, 'disp')")
content.append("ops.recorder('Node', '-file', reaction_file_path, "
"'-node', *supports, '-dof', ctrl_dof, 'reaction')")
# Add lines for saving storey heights
content.append("")
content.append("# Base level coordinate")
content.append(
"base_level = min([ops.nodeCoord(node, 3) for node in supports])"
)
content.append("# Save storey heights")
content.append("with open(storey_heights_file_path, 'w') as file:")
content.append(" for node in floors:")
content.append(
" file.write(f'{ops.nodeCoord(node, 3) - base_level}\\n')"
)
# Add lines for setting the analysis parameters
content.append("")
content.append("# Set some analysis parameters")
content.append(
"max_disp = max_drift * (ops.nodeCoord(ctrl_node, 3) - base_level)"
)
content.append("tol_init = 1.0e-6")
content.append("iter_init = 20")
content.append("ops.wipeAnalysis()")
content.append("ops.system('UmfPack')")
content.append("ops.numberer('RCM')")
content.append("ops.constraints('Transformation')")
content.append("ops.test('EnergyIncr', tol_init, iter_init)")
content.append("ops.integrator('DisplacementControl', ctrl_node, "
"ctrl_dof, dincr)")
content.append("ops.algorithm('Newton', '-initialThenCurrent')")
content.append("ops.analysis('Static')")
# Add lines for performing the analysis
content.append("")
content.append("# Start performing the analysis")
content.append("base_shear = [0]")
content.append("ctrl_disp = [0]")
content.append("ok = 0")
content.append("cont = True")
content.append("while ok == 0 and cont:")
content.append(" # Perform the analysis for a single step with "
"current settings")
content.append(" ok = ops.analyze(1)")
content.append(" if ok != 0: # try other algorithms")
content.append(" ok = _set_algorithm(tol_init, ctrl_node,"
" ctrl_dof, dincr)")
content.append(" if ok != 0: # reduce dincr to an half")
content.append(" ok = _set_algorithm(tol_init, ctrl_node,"
" ctrl_dof, 0.5 * dincr)")
content.append(" if ok != 0: # reduce dincr to a quarter")
content.append(" ok = _set_algorithm(tol_init, ctrl_node,"
" ctrl_dof, 0.25 * dincr)")
content.append(" if ok != 0: # increase tolerance by factor of 10")
content.append(" ok = _set_algorithm(10 * tol_init, ctrl_node,"
" ctrl_dof, 0.25 * dincr)")
content.append(
" if ok != 0: # increase tolerance by factor of 100"
)
content.append(" ok = _set_algorithm(100 * tol_init, ctrl_node,"
" ctrl_dof, 0.25 * dincr)")
content.append("")
content.append(" # Get the base shear force")
content.append(" ops.reactions()")
content.append(" current_disp = ops.nodeDisp(ctrl_node, ctrl_dof)")
content.append(" current_shear = abs(sum([ops.nodeReaction(node, "
"ctrl_dof) for node in supports]))")
content.append(" # Check the base shear limit")
content.append(" within_limits = current_shear >= "
"0.6 * max(base_shear)")
content.append(" # Set continue flag")
content.append(" cont = within_limits and current_disp < max_disp")
content.append(" # Append base shear and control node displacement,"
" including the step reaching max_disp")
content.append(" if ok == 0 and within_limits:")
content.append(" base_shear.append(current_shear)")
content.append(" ctrl_disp.append(current_disp)")
# Add lines for wiping the model
content.append("")
content.append("# Wipe the model")
content.append("ops.wipe()")
# Add lines for return statement
content.append(
"# Return base shear and control node displacement history"
)
content.append("return ctrl_disp, base_shear")
content.append("")
# Add white spaces for method content
content = [' ' + item if item
else item
for item in content]
# Add method definition
method = [
f"def do_nspa_{direction}(max_drift: float = {self.max_drift}, "
f"dincr: float = {self.dincr}) "
"-> tuple[list[float], list[float]]:"
]
content = method + content
return content
def _get_algorithm_py(self) -> List[str]:
"""Generates a list of strings representing the Python function
for setting the algorithm during nspa in OpenSees.
Return
------
List[str]
List of lines representing _set_algorithm method.
"""
content = [
'"""Sets the solution algorithm for NSPA in ops domain.',
"",
"Parameters",
"----------",
"tol : float",
" The tolerance criteria used to check for convergence.",
"ctrl_node : int",
" Tag of control node.",
"ctrl_dof : int",
" Tag of control degrees of freedom, i.e., 1 or 2.",
"dincr : float",
" The displacement increment considered during analysis.",
"iter : int, optional",
" The max number of iterations to check before returning "
"failure.",
" By default 100.",
"",
"Return",
"------",
"int",
" Result of the new analysis step in OpenSees.",
'"""'
]
# Add lines for setting nspa algorithm
content.append("# Set testing and control procedures")
content.append("ops.test('NormDispIncr', tol, iter)")
content.append(
"ops.integrator('DisplacementControl', ctrl_node, ctrl_dof, dincr)"
)
content.append("# Try KrylovNewton")
content.append("ops.algorithm('KrylovNewton')")
content.append("ok = ops.analyze(1)")
content.append("# Try NewtonLineSearch algorithm")
content.append("if ok != 0:")
content.append(" ops.algorithm('NewtonLineSearch', "
"'-InitialInterpolated', 0.8)")
content.append(" ok = ops.analyze(1)")
content.append(
"# Try Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm"
)
content.append("if ok != 0:")
content.append(" ops.algorithm('BFGS')")
content.append(" ok = ops.analyze(1)")
content.append("# Return the analysis result")
# Add lines for setting return statement
content.append("return ok")
content.append("")
# Add white spaces for method content
content = [' ' + item if item
else item
for item in content]
# Add method definition
method = [
"def _set_algorithm(tol: float, ctrl_node: int, ctrl_dof: int,"
" dincr: float, iter: int = 100) -> None:"
]
content = method + content
return content
def _get_nspa_py(self) -> List[str]:
"""Generates a list of strings representing Python functions for
performing nonlinear static pushover (NSPA) analysis in OpenSees.
Return
------
List[str]
List of lines which will be written into nspa.py file.
"""
# Get all nspa methods
algorithm_list = self._get_algorithm_py() + [""]
nspa_x_list = self._get_nspa_dof_py(1) + [""]
nspa_y_list = self._get_nspa_dof_py(2)
# Add imports
imports = ["import openseespy.opensees as ops",
"from pathlib import Path",
"",
"from .model import build_model",
"", ""]
content = imports + algorithm_list + nspa_x_list + nspa_y_list
return content
def _get_run_py(self) -> List[str]:
"""Generates a list of strings representing an example Python script
for running nspa and modal analysis.
Return
------
List[str]
List of lines which will be written into run.py.
"""
content = [
"# Import necessary modules using dynamic import",
"import sys",
"from pathlib import Path",
"import importlib",
"",
"",
"# Determine the package name dynamically",
"current_dir = Path(__file__).resolve().parent",
"parent_dir = current_dir.parent",
"package_name = current_dir.name",
"# Add the parent directory to sys.path",
"sys.path.append(str(parent_dir))",
"# Dynamically import modules from the determined package",
"modal = importlib.import_module(f\"{package_name}.modal\")",
"nspa = importlib.import_module(f\"{package_name}.nspa\")",
"",
"# Perform modal analysis",
"results = modal.do_modal()",
"# Perform nonlinear static pushover analysis in X direction",
"dx, vx = nspa.do_nspa_x()",
"# Perform nonlinear static pushover analysis in Y direction",
"dy, vy = nspa.do_nspa_y()",
""
]
return content
def _get_shared_opstcl(self) -> List[str]:
"""Generate Tcl lines for shared OpenSees domain objects.
Return
------
List[str]
List of lines that define reusable materials, sections, and
geometric transformations that can be referenced by multiple
structural components.
"""
return [
"# Rigid-like material",
f"uniaxialMaterial Elastic {RIGID_MAT} {BIG_VALUE}",
]
def _get_build_tcl(self) -> List[str]:
"""Generates a list of strings representing the Tcl script that
constructs the components of the numerical model in OpenSees.
Return
------
List[str]
List of lines which will be written into model.tcl file.
"""
# Get lines of adding the shared objects
shared = self._get_shared_opstcl()
# Add lines for setting build method
content = [
'# Adds the numerical model to the OpenSees domain',
"",
"# Use the infills flag set by the caller, if any",
"if {![info exists ::infills]} {",
f" set ::infills {self.include_infills}",
"}",
"",
"wipe",
"model BasicBuilder -ndm 3 -ndf 6",
"",
*shared,
"",
"# Define components of foundations",
"source foundations.tcl",
"# Define components of joints",
"source joints.tcl",
"# Define components of floors",
"source floors.tcl",
"# Define components of beams",
"source beams.tcl",
"# Define components of columns",
"source columns.tcl",
"# Perform static analysis under gravity loads",
"source gravity.tcl",
"# Define components of infills (after gravity, unstressed)",
"if {$::infills} {",
" source infills.tcl",
"}",
""
]
return content
def _get_foundations_tcl(self) -> List[str]:
"""Generates a list of strings representing the Tcl script that
constructs the components of the foundation in OpenSees.
Return
------
List[str]
List of lines which will be written into foundations.tcl file.
"""
content = [
'# Add foundation components to ops domain '
'(nodes and constraints)\n'
]
for foundation in self.foundations:
content.extend(foundation.to_tcl())
content.append("")
return content
def _get_floors_tcl(self) -> List[str]:
"""Generates a list of strings representing the Tcl script that
constructs the components of all floors in OpenSees.
Return
------
List[str]
List of lines which will be written into floors.tcl file.
"""
content = [
'# Add floors to ops domain (nodes & diaphrams)\n'
]
for floor in self.floors:
content.extend(floor.to_tcl())
content.append("")
return content
def _get_joints_tcl(self) -> List[str]:
"""Generates a list of strings representing the Tcl script that
constructs the components of all joints in OpenSees.
Return
------
List[str]
List of lines which will be written into joints.tcl file.
"""
content = [
'# Add components of joints to ops domain\n'
]
content.append('# -------------------------------------------------')
content.append('# Add stairs joints to ops domain')
content.append('# -------------------------------------------------\n')
for joint in self.stairs_joints:
content.extend(joint.to_tcl())
content.append("")
content.append('# -------------------------------------------------')
content.append('# Add floor joints to ops domain')
content.append('# -------------------------------------------------\n')
for joint in self.floor_joints:
content.extend(joint.to_tcl())
content.append("")
return content
def _get_beams_tcl(self) -> List[str]:
"""Generates a list of strings representing the Tcl script that
constructs the components of all beams in OpenSees.
Return
------
List[str]
List of lines which will be written into beams.tcl file.
"""
content = [
'# Add components of all beams to ops domain\n'
]
for beam in self.beams:
content.extend(beam.to_tcl())
content.append("")
return content
def _get_columns_tcl(self) -> List[str]:
"""Generates a list of strings representing the Tcl function that
constructs the components of all columns in OpenSees.
Return
------
List[str]
List of lines which will be written into columns.tcl file.
"""
content = [
'# Add components of all columns to ops domain\n'
]
for column in self.columns:
content.extend(column.to_tcl())
content.append("")
return content
def _get_infills_tcl(self) -> List[str]:
"""Generates a list of strings representing the Tcl function that
constructs the components of all the masonry infill walls in OpenSees.
Return
------
List[str]
List of lines which will be written into infills.tcl file.
"""
content = [
'# Add components of all infills to ops domain\n'
]
for infills in self.infills:
content.extend(infills.to_tcl())
content.append("")
return content
def _get_gravity_tcl(self) -> List[str]:
"""Generates a list of strings representing the Tcl script that
constructs the components of the gravity analysis in OpenSees.
Return
------
List[str]
List of lines which will be written into gravity.tcl file.
"""
content = [
"# Perform linear static analysis under gravity loads",
"",
"# Add gravity time-series and load pattern to ops domain",
f"timeSeries Linear {GRAV_TS_TAG}",
f"pattern Plain {GRAV_P_TAG} {GRAV_TS_TAG}" + " {",
""
]
content.append(" # Add beam gravity loads to ops domain")
for beam in self.beams:
content.extend([" " + i for i in beam.to_tcl_grav_loads()])
content.append("")
content.append(" # Add column gravity loads to ops domain")
for column in self.columns:
content.extend([" " + i for i in column.to_tcl_grav_loads()])
content.append("}")
content.append("")
content.append("# Perform gravity analysis and save the model state")
content.append("system UmfPack")
content.append("numberer RCM")
content.append("constraints Transformation")
content.append("test NormDispIncr 1e-08 10 3")
content.append("integrator LoadControl 0.1")
content.append("algorithm Newton")
content.append("analysis Static")
content.append("if {[analyze 10] != 0} {")
content.append(" error \"Gravity analysis did not converge; the "
"model state is invalid for subsequent analyses.\"")
content.append("}")
content.append("loadConst -time 0.0")
content.append("wipeAnalysis")
content.append("")
return content
def _get_modal_tcl(self) -> List[str]:
"""Generates a list of strings representing the Tcl procedure for
modal analysis in OpenSees.
Return
------
List[str]
List of lines which will be written into modal.tcl file.
"""
# Retrained floor nodes
nodes = ' '.join([f"{floor.rnode.tag}" for floor in self.floors])
# Initialize the content list
content = [
"# Perform modal analysis for built OpenSees model",
"#",
"# Parameters",
"# ----------",
"# num_modes : int, optional",
"# Number of modes considered for modal analysis.",
"# By default 3."
"#",
"# Return",
"# ------",
"# dict",
"# Dictionary containing modal properties.",
]
content.append("")
# Add lines for setting the output directory
content.append("# Set output directory")
content.append('set output_directory "Modal-Results"')
content.append("file mkdir $output_directory")
# Add lines for building the numerical model
content.append("# Build the numerical model")
content.append("source model.tcl")
# Add lines for performing eigenvalue analysis
content.append("")
content.append(f"# Perform eigenvalue analysis, with at least "
f"{MIN_EIGEN_MODES} modes (ARPACK is unreliable for a "
"single eigenvalue)")
content.append(f"set eigenVals [eigen [expr {{max($numModes, "
f"{MIN_EIGEN_MODES})}}]]")
content.append("foreach lambda $eigenVals {")
content.append(f" if {{$lambda <= {MIN_EIGEN_VALUE} || "
"$lambda > 1e+300} {")
content.append(' error "Invalid eigenvalues: $eigenVals"')
content.append(" }")
content.append("}")
# Add lines for looping through eigenvectors and saving them
content.append("# Save eigen vectors for retained floor nodes")
content.append(f'set nodes [list {nodes}]')
content.append("# Loop through eigenvectors")
content.append("for {set i 1} {$i <= $numModes} {incr i} {")
content.append(' set modalDisps ""')
content.append(" foreach node $nodes {")
content.append(" # Get the eigenvector as a list of floats")
content.append(" set eigenvector [nodeEigenvector $node $i]")
content.append(" # Join the list into a comma-separated string")
content.append(' set disps [join $eigenvector ", "]')
content.append(" # Append the node tag and the comma-separated"
" eigenvector string")
content.append(' append modalDisps "$node, $disps\\n"')
content.append(" }")
content.append(
' # Write the eigenvector data for the mode to a file')
content.append(
' set report_file '
'[open "$output_directory/EigenVectors_Mode${i}.txt" "w"]'
)
content.append(" puts $report_file $modalDisps")
content.append(" close $report_file")
content.append("}")
# Add lines for saving modal properties
content.append("")
content.append("# Perform modal analysis and save results")
content.append(
'set report_file "$output_directory/ModalProperties.txt"'
)
content.append(
'set results [modalProperties -print -return -file '
'$report_file -unorm]'
)
# Add return statement
content.append("")
content.append("return $results")
# Add Tcl procedure definition
tmp = ["proc do_modal { {numModes 3} } {"]
content = tmp + \
[" " + item if item
else item
for item in content]
# End of procedure
content.append("}")
return content
def _get_nspa_dof_tcl(self, ctrl_dof: Literal[1, 2]) -> List[str]:
"""Generates a list of strings representing the Tcl procedure
for nonlinear static pushover (NSPA) analysis in OpenSees.
Parameters
----------
ctrl_dof : Literal[1, 2]
Control degrees of freedom for loading.
1: X-direction.
2: Y-direction.
Return
------
List[str]
List of lines representing do_nspa_{direction} procedure.
"""
# Get NSPA loading parameters
nodes, loads, ctrl_node = \
self._get_nspa_loading_parameters(ctrl_dof)
# Set support nodes
supports = ' '.join([
f"{found.foundation_node.tag}" for found in self.foundations
])
# Set floor nodes
floors = ' '.join([f"{floor.rnode.tag}" for floor in self.floors])
# Directions per dof
if ctrl_dof == 1:
direction = 'x'
elif ctrl_dof == 2:
direction = 'y'
# Initialize the content list - docstring
content = [
"# Performs nonlinear static pushover analysis (NSPA) "
f"in {direction} direction.",
"#",
"# Parameters",
"# ----------",
"# max_drift : float, optional.",
"# Maximum considered drift value for the control node.",
f"# By default {self.max_drift}",
"# dincr : float, optional.",
"# First displacement increment considered during the "
"analysis.",
f"# By default {self.dincr}.",
"#",
"# Return",
"# ------",
"# ctrl_disp : List[float]",
"# Displacement values of control node.",
"# base_shear : List[float]",
"# Base shear value obtained as sum of the reaction forces.",
""
]
# Add lines for setting the output directory
content.append("# Set output directory")
content.append("set output_directory \"NSPA-Results\"")
content.append("file mkdir $output_directory")
content.append(
'set reaction_file_path "$output_directory'
f'/support_reactions_{direction}.out"')
content.append(
'set disp_file_path "$output_directory'
f'/storey_displacements_{direction}.out"')
content.append(
'set storey_heights_file_path '
'"$output_directory/storey_heights.out"')
# Add lines for building the numerical model
content.append("")
content.append("# Build the numerical model")
content.append("source model.tcl")
# Add lines for setting the time-series and load pattern
content.append("")
content.append("# Add NSPA time-series and load pattern to ops domain")
content.append(f"timeSeries Linear {NSPA_TS_TAG}")
content.append(f"pattern Plain {NSPA_P_TAG} {NSPA_TS_TAG} " + "{")
# Add lines for setting the nspa loads
for node, load_values in zip(nodes, loads):
values_str = ' '.join([f"{val}" for val in load_values])
content.append(f" load {node} {values_str}")
content.append("}")
# Add lines for setting recorders
content.append("")
content.append("# Set recorders")
content.append(f"set ctrl_node {ctrl_node}")
content.append(f"set ctrl_dof {ctrl_dof}")
content.append(f"set supports [list {supports}]")
content.append(f"set floors [list {floors}]")
content.append("recorder Node -file $disp_file_path -node {*}$floors "
"-dof $ctrl_dof disp")
content.append("recorder Node -file $reaction_file_path -node "
"{*}$supports -dof $ctrl_dof reaction")
# Add lines for saving storey heights
content.append("")
content.append("# Base level coordinate")
content.append("set base_level 1.0e12")
content.append("foreach node $supports {")
content.append(" set zCoord [nodeCoord $node 3]")
content.append(" if {$zCoord < $base_level} {")
content.append(" set base_level $zCoord")
content.append(" }")
content.append("}")
content.append("# Save storey heights")
content.append('set file [open $storey_heights_file_path "w"]')
content.append("foreach node $floors {")
content.append(
" puts $file [expr {[nodeCoord $node 3] - $base_level}]"
)
content.append("}")
content.append("close $file")
# Add lines for setting the analysis parameters
content.append("")
content.append("# Set analysis parameters")
content.append(
"set max_disp "
"[expr {$max_drift * ([nodeCoord $ctrl_node 3] - $base_level)}]"
)
content.append("set tol_init 1.0e-6")
content.append("set iter_init 20")
content.append("wipeAnalysis")
content.append("system UmfPack")
content.append("numberer RCM")
content.append("constraints Transformation")
content.append("test EnergyIncr $tol_init $iter_init")
content.append("integrator DisplacementControl $ctrl_node $ctrl_dof "
"$dincr")
content.append("algorithm Newton -initialThenCurrent")
content.append("analysis Static")
# Add lines for performing the analysis
content.append("")
content.append("# Start performing the analysis")
content.append("set max_base_shear 0")
content.append("set ok 0")
content.append("set cont 1")
content.append("set base_shear [list 0.0]")
content.append("set ctrl_disp [list 0.0]")
content.append("while { $ok == 0 && $cont == 1 } {")
content.append(" # Perform the analysis for a single step with "
"current settings")
content.append(" set ok [analyze 1]")
content.append(" # try other algorithms")
content.append(" if { $ok != 0 } {")
content.append(" set ok [_set_algorithm $tol_init $ctrl_node "
"$ctrl_dof $dincr]")
content.append(" }")
content.append(" # reduce dincr to an half")
content.append(" if { $ok != 0 } {")
content.append(" set ok [_set_algorithm $tol_init $ctrl_node "
"$ctrl_dof [expr 0.5 * $dincr]]")
content.append(" }")
content.append(" # reduce dincr to a quarter")
content.append(" if { $ok != 0 } {")
content.append(" set ok [_set_algorithm $tol_init $ctrl_node "
"$ctrl_dof [expr 0.25 * $dincr]]")
content.append(" }")
content.append(" # increase tolerance by factor of 10")
content.append(" if { $ok != 0 } {")
content.append(" set ok [_set_algorithm [expr 10 * $tol_init] "
"$ctrl_node $ctrl_dof [expr 0.25 * $dincr]]")
content.append(" }")
content.append(" # increase tolerance by factor of 100")
content.append(" if { $ok != 0 } {")
content.append(" set ok [_set_algorithm [expr 100 * $tol_init] "
"$ctrl_node $ctrl_dof [expr 0.25 * $dincr]]")
content.append(" }")
content.append("")
content.append(" # Get the base shear force")
content.append(" reactions")
content.append(" set current_disp [nodeDisp $ctrl_node $ctrl_dof]")
content.append(" set sum_reactions 0")
content.append(" foreach foundation_node $supports {")
content.append(
" set reaction [nodeReaction $foundation_node $ctrl_dof]")
content.append(
" set sum_reactions [expr {$sum_reactions + $reaction}]")
content.append(" }")
content.append(" set current_shear [expr {abs($sum_reactions)}]")
content.append(" # Check the base shear limit")
content.append(
" "
"set within_limits [expr {$current_shear >= "
"0.6 * $max_base_shear}]")
content.append(" # Set continue flag")
content.append(
" "
"set cont [expr {$within_limits && "
"($current_disp < $max_disp)}]")
content.append(" # Append base shear and control node displacement,"
" including the step reaching max_disp")
content.append(" if { $ok == 0 && $within_limits } {")
content.append(" lappend base_shear $current_shear")
content.append(" lappend ctrl_disp $current_disp")
content.append(" # Update the maximum recorded shear value")
content.append(
" "
"set max_base_shear [expr {max($max_base_shear, $current_shear)}]"
)
content.append(" }")
content.append("}")
# Add lines for wiping the model
content.append("")
content.append("# Wipe the model")
content.append("wipe")
# Add lines for return statement
content.append(
"# Return base shear and control node displacement history"
)
content.append("return [list $ctrl_disp $base_shear]")
# Add white spaces for method content
content = [' ' + item if item
else item
for item in content]
# Add method definition
var1 = "{ max_drift " + f"{self.max_drift}" + " }"
var2 = "{ dincr " + f"{self.dincr}" + " }"
method = [f"proc do_nspa_{direction}" + " { "
+ f"{var1} {var2}" + " }" + " { "]
content = method + content
# Add method closure bracket
content = content + ["}", ""]
return content
def _get_algorithm_tcl(self) -> List[str]:
"""Generates a list of strings representing the Tcl procedure
for setting the algorithm during nspa in OpenSees.
Return
------
List[str]
List of lines representing Tcl proceudure _set_algorithm.
"""
# Add the procedure content
content = [
"# Sets the solution algorithm for NSPA in OpenSees.",
"#",
"# Parameters",
"# ----------",
"# tol : float",
"# The tolerance criteria used to check for convergence.",
"# ctrl_node : int",
"# Tag of control node.",
"# ctrl_dof : int",
"# Tag of control degrees of freedom, i.e., 1 or 2.",
"# dincr : float",
"# The displacement increment considered during analysis.",
"# iter : float",
"# The max number of iterations to check before returning "
"failure.",
"# By default 100.",
"#",
"# Return",
"# ------",
"# int",
"# Result of the new analysis step in OpenSees.",
"",
""
"# Set testing and control procedures",
"test NormDispIncr $tol $iter",
"integrator DisplacementControl $ctrl_node $ctrl_dof $dincr",
"# Try KrylovNewton",
"algorithm KrylovNewton",
"set ok [analyze 1]",
"# Try NewtonLineSearch algorithm",
"if { $ok != 0 } {",
" algorithm NewtonLineSearch -InitialInterpolated 0.8",
" set ok [analyze 1]",
"}",
"# Try Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm",
"if { $ok != 0 } {",
" algorithm BFGS",
" set ok [analyze 1]",
"}",
"# Return the analysis result",
"return $ok"
]
# Add white spaces for method content
content = [' ' + item if item
else item
for item in content]
# Add method definition
method = ["proc _set_algorithm { tol ctrl_node ctrl_dof dincr "
"{iter 100} } {"]
content = method + content
# Add method closure bracket
content = content + ["}", ""]
return content
def _get_nspa_tcl(self) -> List[str]:
"""Generates a list of strings representing Tcl procedures for
performing nonlinear static pushover (NSPA) analysis in OpenSees.
Return
------
List[str]
List of lines which will be written into nspa.tcl file.
"""
# Get all nspa procedures
algorithm_list = self._get_algorithm_tcl() + [""]
nspa_x_list = self._get_nspa_dof_tcl(1) + [""]
nspa_y_list = self._get_nspa_dof_tcl(2)
content = algorithm_list + nspa_x_list + nspa_y_list
return content
def _get_run_tcl(self) -> List[str]:
"""Generates a list of strings representing an example Python script
for running nspa and modal analysis.
Return
------
List[str]
List of lines which will be written into run.py.
"""
content = [
"source modal.tcl",
"source nspa.tcl",
"",
"# Perform modal analysis",
"set results [do_modal 3]",
"# Perform nonlinear static pushover analysis in X direction",
"lassign [do_nspa_x] dx vx",
"# Perform nonlinear static pushover analysis in Y direction",
"lassign [do_nspa_y] dy vy",
""
]
return content
[docs]
def plot_model(
self, show_nodes: Literal['no', 'yes'] = 'yes', line_width: float = 3,
directory: Optional[str | Path] = None, show: bool = True
) -> None:
"""
Plots the structural model, showing nodes and elements grouped by type
(nodes, beams, columns, infills).
Parameters
----------
show_nodes : Literal['no', 'yes'], optional
A flag to control whether to display the nodes in the plot.
'yes' to show the nodes, 'no' to hide them. By default 'yes'.
line_width : float, optional
Specifies the line width used to draw the elements in the plot.
By default 3.
directory : str | Path | None, optional
Directory to save an image of the model. If None, the image will
not be saved. By default None.
show : bool, optional
Flag for showing the figure in an interactive window,
by default True.
"""
# Set the group elements
beams = [beam.design.line.tag for beam in self.beams]
columns = [column.design.line.tag for column in self.columns]
infills = []
if self.include_infills:
for infill in self.infills:
s1 = sub_tag(infill.design.rectangle.tag, 1)
s2 = sub_tag(infill.design.rectangle.tag, 2)
infills.append(s1)
infills.append(s2)
groups = [
[beams, columns, infills],
["red", "blue", "green"]
]
# Build the model
self.build()
# Path to the file (without the file extension)
if directory:
filename = str(Path(directory) / 'model_view.html')
if not Path.exists(Path(directory)):
make_dir(directory)
else:
filename = None
# Plot the model
pl.plot_model(show_nodes=show_nodes, ele_groups=groups, show=show,
line_width=line_width, filename=filename)
[docs]
def plot_mode_shape(
self, mode_number: int = 1, scale: float = 100, line_width: float = 3,
contour: Optional[Literal['x', 'y', 'z']] = None,
set_view: Literal['xy', 'yz', 'xz', '3D'] = '3D',
directory: Optional[str | Path] = None, show: bool = True
) -> None:
"""
Plots the mode shape of the structure for a given mode number,
scaled for visual clarity.
Parameters
----------
mode_number : int, optional
Specifies the mode number to be plotted.
By default 1.
scale : float, optional
A scaling factor to exaggerate the mode shape for better
visualization. By default 100.
line_width : float, optional
Specifies the line width used to draw the mode shape.
By default 3.
contour : Literal['x', 'y', 'z'] | None, optional
Contours of displacement in x, y, or z.
By default None.
set_view : Literal['xy', 'yz', 'xz', '3D'], optional
Camera view plane. By default '3D' (isometric).
directory : str | Path | None, optional
Directory to save an image of the model. If None, the image will
not be saved. By default None.
show : bool, optional
Flag for showing the figure in an interactive window,
by default True.
"""
# Build the model
self.build()
# Path to the file
if directory:
filename = str(Path(directory) / f'mode_{mode_number}_shape.html')
if not Path.exists(Path(directory)):
make_dir(directory)
else:
filename = None
# Plot the mode shape
pl.plot_mode_shape(
mode_number=mode_number, scale=scale, line_width=line_width,
contour=contour, set_view=set_view, filename=filename, show=show
)