"""
Define the functions for defining and manipulating a model.
"""
from math import sqrt, pi
from numpy import array, zeros, dot, mean, concatenate, float64, int32, inf
import scipy
from scipy.linalg import solve, eigh
from collections import namedtuple
from pystran import truss, beam, spring, rigid
from numbers import Integral
[docs]
def create(dim=2):
"""
Create a new model.
Parameters
----------
dim
Supply the dimension of the model (2 or 3).
Returns
-------
dict
The model is represented as a dictionary. The keys of
the dictionary are described in the documentation of the functions that manipulate the model.
Initially, the model contains only values for the the key ``dim`` that gives the dimension of the model,
and for the key ``freedoms``, which is an enumeration of the degrees of freedom at a joint.
2D models have three degrees of freedom at a joint: two translations (`m['freedoms'].U1`,
`m['freedoms'].U2`) and one rotation (`m['freedoms'].UR3`).
3D models have six degrees of freedom at a joint: three translations (`m['freedoms'].U1`,
`m['freedoms'].U2`, `m['freedoms'].U3`) and three rotations (`m['freedoms'].UR1`,
`m['freedoms'].UR2`, `m['freedoms'].UR3`).
See Also
--------
:func:`add_joint`
:func:`add_truss_member`
:func:`add_beam_member`
:func:`add_rigid_link_member`
:func:`add_spring_member`
:func:`add_support`
:func:`add_load`
:func:`add_mass`
:func:`add_dof_links`
"""
m = {}
m["dim"] = dim # Dimension of the model
if m["dim"] == 2:
Freedoms = namedtuple("Freedoms", ["U1", "U2", "UR3",
"TRANSLATION_DOFS",
"ROTATION_DOFS",
"ALL_DOFS"])
m["freedoms"] = Freedoms(U1=0, U2=1, UR3=2,
TRANSLATION_DOFS=(0, 1),
ROTATION_DOFS=(2,),
ALL_DOFS=(0, 1, 2))
else:
Freedoms = namedtuple("Freedoms", ["U1", "U2", "U3", "UR1", "UR2", "UR3",
"TRANSLATION_DOFS",
"ROTATION_DOFS",
"ALL_DOFS"])
m["freedoms"] = Freedoms(U1=0, U2=1, U3=2, UR1=3, UR2=4, UR3=5,
TRANSLATION_DOFS=(0, 1, 2),
ROTATION_DOFS=(3, 4, 5),
ALL_DOFS=(0, 1, 2, 3, 4, 5))
return m
[docs]
def add_joint(m, jid, coordinates, dof=None):
"""
Add a joint to the model.
Parameters
----------
m
Model.
jid
The joint identifier, which must be unique, but can be anything
that is a legal dictionary key (integer, string, ...),
coordinates
The list (or a tuple) of coordinates of the joint; the
input is converted to an array.
dof
Optional: the degrees of freedom of the joint as a list (or a
tuple). If provided, do not use :func:`number_dofs` later on.
Returns
-------
dict
Newly created joint.
"""
if "joints" not in m:
m["joints"] = {}
if jid in m["joints"]:
raise RuntimeError("Joint already exists")
coordinates = array(coordinates, dtype=float64)
if coordinates.shape != (m["dim"],):
raise RuntimeError("Coordinate dimension mismatch")
m["joints"][jid] = {"jid": jid, "coordinates": coordinates}
if dof is not None:
m["joints"][jid]["dof"] = array(dof, dtype=int32)
return m["joints"][jid]
[docs]
def add_truss_member(m, mid, connectivity, sect):
"""
Add a truss member to the model.
Parameters
----------
m
Model.
mid
The member identifier, which must be unique, but can be anything
that is a legal dictionary key (integer, string, ...).
connectivity
The list (or a tuple) of the joint identifiers.
sect
Section of appropriate type (i.e. a truss section).
Returns
-------
dict
Newly created member.
"""
if "truss_members" not in m:
m["truss_members"] = {}
if mid in m["truss_members"]:
raise RuntimeError("Truss member already exists")
m["truss_members"][mid] = {
"mid": mid,
"connectivity": connectivity,
"section": sect,
}
return m["truss_members"][mid]
[docs]
def add_beam_member(m, mid, connectivity, sect):
"""
Add a beam member to the model.
Parameters
----------
m
Model.
mid
The member identifier, which must be unique, but can be anything
that is a legal dictionary key (integer, string, ...).
connectivity
The list (or a tuple) of the joint identifiers.
sect
Section of appropriate type (2d or 3d beam section).
Returns
-------
dict
Newly created member.
See Also
--------
:func:`section.beam_2d_section`
:func:`section.beam_3d_section`
"""
if "beam_members" not in m:
m["beam_members"] = {}
if mid in m["beam_members"]:
raise RuntimeError("Beam member already exists")
m["beam_members"][mid] = {
"mid": mid,
"connectivity": connectivity,
"section": sect,
}
return m["beam_members"][mid]
[docs]
def add_rigid_link_member(m, mid, connectivity, sect):
"""
Add a rigid link member to the model.
Parameters
----------
m
Model.
mid
The member identifier, which must be unique, but can be anything
that is a legal dictionary key (integer, string, ...),
connectivity
The list (or a tuple) of the joint identifiers; the
first is the *master* (its motion determines the motion of the
subordinate), the second is the *subordinate* (its motion follows that of
the master).
sect
Section of appropriate type (i.e. rigid link section).
Returns
-------
dict
Newly created member.
See Also
--------
:func:`rigid_link_section`
"""
if "rigid_link_members" not in m:
m["rigid_link_members"] = {}
if mid in m["rigid_link_members"]:
raise RuntimeError("Rigid link member already exists")
m["rigid_link_members"][mid] = {
"mid": mid,
"connectivity": connectivity,
"section": sect,
}
return m["rigid_link_members"][mid]
[docs]
def add_spring_member(m, mid, connectivity, sect):
"""
Add a spring member to the model.
Parameters
----------
m
Model.
mid
The member identifier, which must be unique, but can be anything
that is a legal dictionary key (integer, string, ...),
connectivity
The list (or a tuple) of the joint identifiers.
sect
Section of appropriate type (i.e. a spring section).
Returns
-------
dict
Newly created member.
See Also
--------
:func:`spring_section`
"""
if "spring_members" not in m:
m["spring_members"] = {}
if mid in m["spring_members"]:
raise RuntimeError("Spring member already exists")
if sect['kind'] == "extension":
dof = m["freedoms"].TRANSLATION_DOFS
else: # torsion
dof = m["freedoms"].ROTATION_DOFS
m["spring_members"][mid] = {
"mid": mid,
"connectivity": connectivity,
"section": sect,
"dofkind": dof
}
return m["spring_members"][mid]
def _dof_is_int(dof):
return isinstance(dof, Integral)
[docs]
def add_support(j, dof, value=0.0):
"""
Add a support to a joint.
Parameters
----------
j
The joint (obtained from the model as ``m["joints"][jid]``).
dof
The degree of freedom (0, 1, ...). Refer to the model key ``'freedoms'``.
value
The signed magnitude of the support motion (default is zero).
Returns
-------
None
"""
if "supports" not in j:
j["supports"] = {}
dim = len(j["coordinates"])
if not _dof_is_int(dof):
for d in dof:
j["supports"][d] = value
else:
j["supports"][dof] = value
[docs]
def add_load(j, dof, value):
"""
Add a load to a joint.
Parameters
----------
j
The joint (obtained from the model as ``m["joints"][jid]``).
dof
The degree of freedom (0, 1, ...). Refer to the model key ``'freedoms'``.
value
The signed magnitude of the load.
Returns
-------
None
"""
if "loads" not in j:
j["loads"] = {}
if dof not in j["loads"]:
j["loads"][dof] = 0.0
j["loads"][dof] += value
[docs]
def add_mass(j, dof, value):
"""
Add a mass to a joint.
Parameters
----------
j
The joint (obtained from the model as ``m["joints"][jid]``).
dof
The degree of freedom (0, 1, ...). Refer to the model key ``'freedoms'``.
value
The magnitude of the added mass.
Returns
-------
None
"""
if "masses" not in j:
j["masses"] = {}
j["masses"][dof] = value
[docs]
def add_dof_links(m, jids, dof):
"""
Add degree-of-freedom links between all joints in the list ``jids`` in the
direction ``dof``.
Parameters
----------
m
The model.
jids
The list of joint identifiers.
dof
The degree of freedom at which the joints are to be linked.
Returns
-------
None
"""
# Now add the mutual links between the joints
if 'joints' not in m:
raise RuntimeError("No joints in the model")
for jid1 in jids:
for jid2 in jids:
if jid1 != jid2:
j1 = m["joints"][jid1]
if "links" not in j1:
j1["links"] = {}
if jid2 not in j1["links"]:
j1["links"][jid2] = []
if _dof_is_int(dof):
j1["links"][jid2].append(dof)
else:
for d in dof:
j1["links"][jid2].append(d)
[docs]
def bounding_box(m):
"""
Compute the bounding box of the model.
Parameters
----------
m
The model.
Returns
-------
array
Array of the lower and upper ranges (i.e. the bounding box that
encloses all the joints).
"""
dim = m["dim"]
box = array(
concatenate([[inf for i in range(dim)], [-inf for i in range(dim)]])
)
for j in m["joints"].values():
cj = j["coordinates"]
for i, v in enumerate(cj):
box[i] = min(box[i], v)
box[i + dim] = max(box[i + dim], v)
return box
[docs]
def characteristic_dimension(m):
"""
Compute the characteristic dimension of the model.
This is the average of the dimensions of the bounding box.
Parameters
----------
m
The model.
Returns
-------
float
Characteristic dimension.
"""
dim = m["dim"]
box = bounding_box(m)
dl = [box[i + dim] - box[i] for i in range(dim)]
return mean(array(dl))
def _copy_dof_num_to_linked(m, j, d, n):
if "links" in j:
for k in j["links"].keys():
o = m["joints"][k]
if d in o["links"][j["jid"]]:
o["dof"][d] = n
def _have_rotations(m):
with_rotations = "beam_members" in m and m["beam_members"]
if with_rotations:
return True
f = m['freedoms']
for j in m["joints"].values():
if "supports" in j and j["supports"]:
for dof in j["supports"].keys():
if (m['dim'] == 2) and dof == f.UR3:
return True
if (m['dim'] == 3) and (dof == f.UR1 or dof == f.UR2 or dof == f.UR3 or dof == f.UR3):
return True
return False
[docs]
def ndof_per_joint(m):
"""
How many degrees of freedom are there per joint?
Parameters
----------
m
The model.
Returns
-------
int
How many degrees of freedom are there per joint? Depends on the space
dimension of the model and the presence or absence of beams.
"""
ndpn = m["dim"]
with_rotations = _have_rotations(m)
if with_rotations:
if m["dim"] == 2:
ndpn = 3
else:
ndpn = 6
return ndpn
[docs]
def number_dofs(m):
"""
Number degrees of freedom.
All current information about degrees of freedom will be replaced when this
function is done.
After this function returns, ``m["nfreedof"]`` will be the number of free
degrees of freedom, and ``m["ntotaldof"]`` will be the total number of degrees
of freedom.
The degrees of freedom are numbered in the order of free and then
prescribed.
Parameters
----------
m
The model.
Returns
-------
None
"""
if "joints" not in m:
raise RuntimeError("No joints in the model")
# Determine the number of degrees of freedom per joint
ndpn = ndof_per_joint(m)
# Generate arrays for storing the degrees of freedom
for j in m["joints"].values():
if "dof" not in j:
j["dof"] = zeros((ndpn,), dtype=int32)
j["dof"][:] = -1 # -1 means not yet numbered
# For each linked pair of joints, make sure they share the same supports
for j in m["joints"].values():
if "links" in j and "supports" in j:
for k in j["links"].keys():
o = m["joints"][k]
if not "supports" in o:
o["supports"] = j["supports"].copy()
if o["supports"] != j["supports"]:
raise RuntimeError("Linked joints must have the same supports")
# Number the free degrees of freedom first
n = 0
for j in m["joints"].values():
for d in range(ndpn):
if ("supports" not in j) or (d not in j["supports"]):
if j["dof"][d] < 0:
j["dof"][d] = n
_copy_dof_num_to_linked(m, j, d, n)
n += 1
m["nfreedof"] = n
# Number all prescribed degrees of freedom
for j in m["joints"].values():
for d in range(ndpn):
if "supports" in j and d in j["supports"]:
if j["dof"][d] < 0:
j["dof"][d] = n
_copy_dof_num_to_linked(m, j, d, n)
n += 1
m["ntotaldof"] = n
def _build_stiffness_matrix(m):
nt = m["ntotaldof"]
# Assemble global stiffness matrix and mass matrix
K = zeros((nt, nt))
if "truss_members" in m:
for member in m["truss_members"].values():
connectivity = member["connectivity"]
i, j = m["joints"][connectivity[0]], m["joints"][connectivity[1]]
truss.assemble_stiffness(K, member, i, j)
if "beam_members" in m:
for member in m["beam_members"].values():
connectivity = member["connectivity"]
i, j = m["joints"][connectivity[0]], m["joints"][connectivity[1]]
beam.assemble_stiffness(K, member, i, j)
if "rigid_link_members" in m:
for member in m["rigid_link_members"].values():
connectivity = member["connectivity"]
i, j = m["joints"][connectivity[0]], m["joints"][connectivity[1]]
rigid.assemble_stiffness(K, member, i, j)
if "spring_members" in m:
for member in m["spring_members"].values():
connectivity = member["connectivity"]
i, j = m["joints"][connectivity[0]], m["joints"][connectivity[1]]
spring.assemble_stiffness(K, member, i, j)
return K
def _build_mass_matrix(m):
nt = m["ntotaldof"]
M = zeros((nt, nt))
if "truss_members" in m:
for member in m["truss_members"].values():
connectivity = member["connectivity"]
i, j = m["joints"][connectivity[0]], m["joints"][connectivity[1]]
truss.assemble_mass(M, member, i, j)
if "beam_members" in m:
for member in m["beam_members"].values():
connectivity = member["connectivity"]
i, j = m["joints"][connectivity[0]], m["joints"][connectivity[1]]
beam.assemble_mass(M, member, i, j)
for j in m["joints"].values():
if "masses" in j:
for dof, value in j["masses"].items():
if _dof_is_int(dof):
gr = j["dof"][dof]
M[gr, gr] += value
else:
for d in dof:
gr = j["dof"][d]
M[gr, gr] += value
return M
[docs]
def solve_statics(m):
r"""
Solve the static equilibrium of the discrete model.
This function solves the equation of static equilibrium
.. math::
K \cdot U = F
Here :math:`K` is the stiffness matrix, :math:`U` is the displacement
vector, and :math:`F` is the vector of forces acting on the joints.
Note that the degrees of freedom can be partitioned into 'free' (unknown)
and 'data' (given, i.e. prescribed).
.. math::
\left[ \begin{array}{cc}
K_{ff} & K_{fd} \\
K_{df} & K_{dd} \\
\end{array}\right] \cdot \left[ \begin{array}{cc}
U_{f} \\
U_{d} \\
\end{array}\right] = \left[ \begin{array}{cc}
L_{f} \\
L_{d} + R\\
\end{array}\right]
Here :math:`L_f` is the vector of active loads applied to the free degrees
of freedom, and :math:`L_d` is the vector of active loads applied to the
data degrees of freedom. The reactions :math:`R` due to supports act on the
prescribed (data) degrees of freedom.
The system of equations is solved for the free degrees of freedom as
.. math::
K_{ff} \cdot U_{f} = -K_{fd} \cdot U_{d} +L_{f}
Note: :func:`number_dofs` must be called before this function to number the
degrees of freedom, automatically. Alternatively, the user may specify the
numbers of the degrees of freedom when defining the joints: the manual way.
Parameters
----------
m
The model.
Returns
-------
None
See Also
--------
:func:`number_dofs`
"""
if not ("ntotaldof" in m) or m["ntotaldof"] <= 0:
raise RuntimeError(
"No degrees of freedom: the numbers of degrees of freedom need to be generated"
)
if not ("nfreedof" in m) or m["nfreedof"] <= 0:
raise RuntimeError("No free degrees of freedom: nothing to compute")
nt, nf = m["ntotaldof"], m["nfreedof"]
# Assemble global stiffness matrix
K = _build_stiffness_matrix(m)
m["K"] = K
# Compute the active load vector
F = zeros(nt)
for joint in m["joints"].values():
if "loads" in joint:
for dof, value in joint["loads"].items():
gr = joint["dof"][dof]
F[gr] += value
m["F"] = F
# Set the prescribed displacements in the displacement vector.
U = zeros(m["ntotaldof"])
for joint in m["joints"].values():
if "supports" in joint:
for dof, value in joint["supports"].items():
if value != 0.0:
gr = joint["dof"][dof]
U[gr] = value
# # Solve for displacements
U[0:nf] = solve(K[0:nf, 0:nf], F[0:nf] - dot(K[0:nf, nf:nt], U[nf:nt]))
m["U"] = U
# # Assign displacements back to joints
for joint in m["joints"].values():
joint["displacements"] = U[joint["dof"]]
[docs]
def statics_reactions(m):
r"""
Compute the reactions in the static equilibrium of the discrete model.
The partitioned system of the balance equations reads
.. math::
\left[ \begin{array}{cc}
K_{ff} & K_{fd} \\
K_{df} & K_{dd} \\
\end{array}\right] \cdot \left[ \begin{array}{cc}
U_{f} \\
U_{d} \\
\end{array}\right] = \left[ \begin{array}{cc}
L_{f} \\
L_{d} + R\\
\end{array}\right]
Here :math:`L_f` is the vector of active loads applied to the free degrees
of freedom, and :math:`L_d` is the vector of active loads applied to the
data degrees of freedom. The reactions :math:`R` due to supports act on the
prescribed (data) degrees of freedom.
The system of equations is solved for the reactions as
.. math::
R = K_{ff} \cdot U_{f} + K_{fd} \cdot U_{d} -L_{d}
once :math:`U_f` has been solved for in the :func:`solve_statics` step.
The reactions are distributed to the joints, and can be retrieved from
individual joint dictionaries ``j`` as ``j['reactions']``.
Parameters
----------
m
The model.
Returns
-------
None
See Also
--------
:func:`solve_statics`
"""
if not ("K" in m):
raise RuntimeError(
"No stiffness matrix: the stiffness matrix needs to be generated by calling solve_statics"
)
K = m["K"]
U = m["U"]
F = m["F"]
# Compute reactions from the partitioned stiffness matrix and the
# partitioned displacement vector
# R = dot(K[nf:nt, 0:nf], U[0:nf]) + dot(K[nf:nt, nf:nt], U[nf:nt]) - F[nf:nt]
# For convenience when working
# with degrees of freedom, we compute this product and only use the rows
# corresponding to fixed the degrees of freedom.
R = dot(K, U) - F
for joint in m["joints"].values():
if "supports" in joint:
reactions = {}
for dof, _ in joint["supports"].items():
gr = joint["dof"][dof]
reactions[dof] = R[gr]
joint["reactions"] = reactions
[docs]
def solve_free_vibration(m, freqshift=0.0):
r"""
Solve the free vibration of the discrete model.
The free vibration eigenvalue problem is solved for the eigenvalues and
eigenvectors (can be retrieved as ``m["eigvals"]`` and ``m["eigvecs"]``). The
frequencies are computed from the eigenvalues (can be retrieved as
``m["frequencies"]``).
The equation of free vibration is
.. math::
K \cdot V = \omega^2 M \cdot V
where :math:`M` is the mass matrix, :math:`V` is the eigenvector, and
:math:`\omega` is the angular frequency.
:func:`number_dofs` must be called before this function.
Parameters
----------
m
The model.
freqshift
Optional: a frequency shift to apply to the eigenvalues, in order to
compute the frequencies around a certain value. The shifted eigenvalue
problem is (with :math:`\bar\omega=2\pi\text{freqshift}` being the shifted angular frequency)
.. math::
(K + \bar\omega^2 M) \cdot V = (\omega^2 - \bar\omega^2) M \cdot V
Returns
-------
None
See Also
--------
:func:`number_dofs`
"""
if not ("ntotaldof" in m) or m["ntotaldof"] <= 0:
raise RuntimeError(
"No degrees of freedom: the numbers of degrees of freedom need to be generated"
)
if not ("nfreedof" in m) or m["nfreedof"] <= 0:
raise RuntimeError("No free degrees of freedom: nothing to compute")
nt, nf = m["ntotaldof"], m["nfreedof"]
# Assemble global stiffness matrix and mass matrix
K = _build_stiffness_matrix(m)
M = _build_mass_matrix(m)
m["K"] = K
m["M"] = M
U = zeros(m["ntotaldof"])
for joint in m["joints"].values():
if "supports" in joint:
for dof, _ in joint["supports"].items():
gr = joint["dof"][dof]
U[gr] = 0.0
m["U"] = U
# Solve the eigenvalue problem. Potentially with shifting for better convergence around a certain frequency.
Kff = K[0:nf, 0:nf]
Mff = M[0:nf, 0:nf]
if freqshift != 0.0:
baromega = (2 * pi * freqshift)
eigvals, eigvecs = eigh(Kff + baromega**2 * Mff, Mff)
eigvals -= baromega**2
else:
eigvals, eigvecs = eigh(Kff, Mff)
m["eigvals"] = eigvals
m["frequencies"] = [sqrt(abs(ev)) / 2 / pi for ev in eigvals]
m["eigvecs"] = eigvecs
[docs]
def set_solution(m, V):
"""
Set the displacement solution from a vector.
Parameters
----------
m
The model.
V
The displacement vector. Either of length ``m["nfreedof"]`` for only the
free degrees of freedom, or of length ``m["ntotaldof"]`` for the total
number of degrees of freedom.
Returns
-------
None
See Also
--------
:func:`number_dofs`
"""
if not ("ntotaldof" in m) or m["ntotaldof"] <= 0:
raise RuntimeError(
"No degrees of freedom: the numbers of degrees of freedom need to be generated"
)
if not ("nfreedof" in m) or m["nfreedof"] <= 0:
raise RuntimeError("No free degrees of freedom: nothing to compute")
nt, nf = m["ntotaldof"], m["nfreedof"]
if len(V) == nf:
m["U"][0:nf] = V
elif len(V) == nt:
m["U"][0:nt] = V
else:
raise RuntimeError("Invalid vector length")
for joint in m["joints"].values():
joint["displacements"] = m["U"][joint["dof"]]
[docs]
def free_body_check(m):
"""
Check the balance of the structure as a free body.
All the active forces and moments together with the reactions at all the
supports should sum to zero.
:func:`statics_reactions` must be called before this function as this calculation
relies on the presence of reactions at the joints.
Parameters
----------
m
The model.
Returns
-------
array
Array of resultant forces and moments.
See Also
--------
:func:`statics_reactions`
"""
if not ("ntotaldof" in m) or m["ntotaldof"] <= 0:
raise RuntimeError(
"No degrees of freedom: the numbers of degrees of freedom need to be generated"
)
if not ("nfreedof" in m) or m["nfreedof"] <= 0:
raise RuntimeError("No free degrees of freedom: nothing to compute")
nt, nf = m["ntotaldof"], m["nfreedof"]
if not ("K" in m):
raise RuntimeError(
"No stiffness matrix: the stiffness matrix needs to be generated by calling solve_statics"
)
if m["dim"] == 2:
nrbm = 3 # Number of rigid body modes: assume 2 translations, 1 rotation
MZ = 2
allforces = zeros(nrbm)
for joint in m["joints"].values():
c = joint["coordinates"]
x, y = c[0], c[1]
if "loads" in joint:
for dof, value in joint["loads"].items():
if dof < MZ:
# Add contributions of forces to the moment
if dof == 0:
allforces[MZ] += -value * y
elif dof == 1:
allforces[MZ] += +value * x
else:
# Add contributions of forces and moments
allforces[dof] += value
if "reactions" in joint:
for dof, value in joint["reactions"].items():
if dof < MZ:
# Add contributions of forces to the moment
if dof == 0:
allforces[MZ] += -value * y
elif dof == 1:
allforces[MZ] += +value * x
else:
# Add contributions of forces and moments
allforces[dof] += value
return allforces
else:
nrbm = 6 # Number of rigid body modes: assume 3 translations, 3 rotations
MX, MY, MZ = 3, 4, 5
allforces = zeros(nrbm)
for joint in m["joints"].values():
c = joint["coordinates"]
x, y, z = c[0], c[1], c[2]
if "loads" in joint:
for dof, value in joint["loads"].items():
if dof < MX:
# Add contributions of forces to the moment
if dof == 0:
allforces[MY] += +value * z
allforces[MZ] += -value * y
elif dof == 1:
allforces[MX] += -value * z
allforces[MZ] += +value * x
else:
allforces[MY] += -value * x
allforces[MX] += +value * y
else:
# Add contributions of forces and moments
allforces[dof] += value
if "reactions" in joint:
for dof, value in joint["reactions"].items():
if dof < MX:
# Add contributions of forces to the moment
if dof == 0:
allforces[MY] += +value * z
allforces[MZ] += -value * y
elif dof == 1:
allforces[MX] += -value * z
allforces[MZ] += +value * x
else:
allforces[MY] += -value * x
allforces[MX] += +value * y
else:
# Add contributions of forces and moments
allforces[dof] += value
return allforces
[docs]
def refine_member(m, mid, n):
"""
Refine a beam member by replacing it with ``n`` new members.
The new joints are numbered starting from zero, and the joint identifier is
composed of the member identifier plus the serial number of the new joint.
The new member identifiers are stored under the key ``"descendants"`` in
the refined member. The refined member is removed from the list of beam
members.
Parameters
----------
m
The model.
mid
The identifier of the member to be refined.
n
The number of new beam members to replace the old member with.
Returns
-------
None
"""
if n < 2:
raise RuntimeError("Number of new members must be at least 2")
if not ("beam_members" in m) or mid not in m["beam_members"]:
raise RuntimeError("Beam member to be refined does not exist")
member = m["beam_members"][mid]
connectivity = member["connectivity"]
i, j = m["joints"][connectivity[0]], m["joints"][connectivity[1]]
ci, cj = i["coordinates"], j["coordinates"]
# Store the descendants
descendants = []
# First replacement member
start = -1.0 + 2.0 / n
c = (-1 + start) / (-2) * ci + (1 + start) / 2 * cj
newjid = str(mid) + "j" + "0"
add_joint(m, newjid, c)
newmid = str(mid) + "m" + "0"
add_beam_member(m, newmid, [i["jid"], newjid], member["section"])
descendants.append(newmid)
prevjid = newjid
for k in range(n - 2):
start += 2.0 / n
c = (-1 + start) / (-2) * ci + (1 + start) / 2 * cj
newjid = str(mid) + "j" + str(k + 1)
add_joint(m, newjid, c)
newmid = str(mid) + "m" + str(k + 1)
add_beam_member(m, newmid, [prevjid, newjid], member["section"])
descendants.append(newmid)
prevjid = newjid
# Last replacement member
newmid = str(mid) + "m" + str(n - 1)
add_beam_member(m, newmid, [newjid, j["jid"]], member["section"])
descendants.append(newmid)
# Remember the provenance of the new members
member["descendants"] = descendants
# Remove the old member
del m["beam_members"][mid]
[docs]
def remove_loads(m):
"""
Remove all the nodal loads in the model.
Parameters
----------
m
The model.
"""
for joint in m["joints"].values():
if "loads" in joint:
joint["loads"] = {}
return None
[docs]
def remove_supports(m):
"""
Remove all the nodal supports in the model.
Parameters
----------
m
The model.
"""
for joint in m["joints"].values():
if "supports" in joint:
joint["supports"] = {}
return None
[docs]
def volume(m):
"""
Compute the total volume of the members in the model.
Parameters
----------
m
The model.
Returns
-------
float
Total volume of the members in the model.
"""
total_volume = 0.0
if "truss_members" in m:
for member in m["truss_members"].values():
connectivity = member["connectivity"]
i, j = m["joints"][connectivity[0]], m["joints"][connectivity[1]]
total_volume += truss.truss_volume(member, i, j)
if "beam_members" in m:
for member in m["beam_members"].values():
connectivity = member["connectivity"]
i, j = m["joints"][connectivity[0]], m["joints"][connectivity[1]]
total_volume += beam.beam_volume(member, i, j)
return total_volume