aframe
aframe: a differentiable linear 3D beam/frame solver written in CSDL.
aframe is a linear finite-element solver for 3D beams and frames, written in CSDL. The code is written entirely in CSDL, so it is fully differentiable and well-suited for gradient-based design optimization and analysis.
📖 Documentation: lsdolab.github.io/aframe
Features
- 3D Euler-Bernoulli beam elements with 6 dofs per node and consistent mass matrices
- Frames of many beams connected by rigid joints, with fixed and pinned supports
- Nodal forces and moments, and inertial loads from a frame acceleration (e.g. gravity) and point masses
- Cross-sections: tube, box (wing box), solid circle, ellipse, or your own properties
- Stress recovery (von Mises) for tube and box sections
- Exact derivatives through CSDL, with the NumPy (inline) or JAX backend
Installation
aframe requires Python 3.9+ and CSDL_alpha, which is not on PyPI, so install it first:
pip install git+https://github.com/LSDOlab/CSDL_alpha.git
git clone https://github.com/LSDOlab/aframe.git
pip install -e ./aframe
Add extras as needed: pip install -e "./aframe[plot]" for the PyVista/Matplotlib plotting
helpers, [dev] for the tests and [docs] to build the documentation. Install jax to use
CSDL's JAX backend.
Quick start: a cantilever beam
A 10 m aluminum tube, clamped at the root, with a 1 kN tip load:
import numpy as np
import csdl_alpha as csdl
import aframe as af
# record the model; inline=True evaluates values as the model is built
recorder = csdl.Recorder(inline=True)
recorder.start()
# a 10 m cantilever along x with 21 nodes (20 elements)
num_nodes = 21
mesh = np.zeros((num_nodes, 3))
mesh[:, 0] = np.linspace(0, 10, num_nodes)
# aluminum tube cross-section, one value per element
radius = csdl.Variable(value=np.full(num_nodes - 1, 0.1), name='radius')
thickness = csdl.Variable(value=np.full(num_nodes - 1, 0.005), name='thickness')
cs = af.CSTube(radius=radius, thickness=thickness)
beam = af.Beam(name='cantilever', mesh=csdl.Variable(value=mesh),
E=69e9, G=26e9, density=2700, cs=cs)
beam.fix(0) # clamp the root node
# nodal loads [Fx, Fy, Fz, Mx, My, Mz]: 1 kN downward at the tip
loads = np.zeros((num_nodes, 6))
loads[-1, 2] = -1000
beam.add_load(csdl.Variable(value=loads))
frame = af.Frame(beams=[beam])
frame.solve()
tip_deflection = frame.displacement['cantilever'][-1, 2]
stress = frame.compute_stress()['cantilever'] # von Mises, one value per element
# derivatives of any output with respect to any input
d_deflection_d_radius = csdl.derivative(tip_deflection, radius)
recorder.stop()
print('tip deflection (m):', tip_deflection.value)
print('max stress (Pa):', stress.value.max())
tip deflection (m): [-0.33159692]
max stress (Pa): 66924549.00143065
Both match beam theory: the tip deflection is PL³/3EI, and the maximum stress is the
bending stress M r / I at the midpoint of the root element (stresses are evaluated at
element midpoints). d_deflection_d_radius has shape (1, 20): the sensitivity of the
tip deflection to the radius of each element.
Building models
Beams
af.Beam(name, mesh, E, G, density, cs, z=False) is a beam discretized into
num_nodes - 1 elements between the nodes of mesh (shape (num_nodes, 3)). E, G
and density are the Young's modulus, shear modulus and density; use any consistent
units (the examples use SI). Beam names must be unique within a frame, since results
are keyed by name.
Each element's local x axis points along the element. The standard construction of the
local y and z axes is singular for elements parallel to the global z axis, so vertical
beams need z=True.
Cross-sections
Section properties are per element, i.e. arrays of shape (num_nodes - 1,), and can be
design variables.
| Cross-section | Parameters | stress() |
|---|---|---|
af.CSTube |
radius, thickness |
von Mises at the outer fiber, shape (num_elements,) |
af.CSBox |
ttop, tbot, tweb, height, width |
von Mises at 4 corners and the web mid-height, shape (num_elements, 5) |
af.CSBox.from_properties |
box dimensions plus area, ix, iy, iz |
as CSBox |
af.CSCircle |
radius |
none |
af.CSEllipse |
semi_major_axis, semi_minor_axis |
none |
af.CrossSection |
area, ix, iy, iz |
none |
ix is the torsion constant and iy, iz are the second moments of area about the
element's local y and z axes. To add a section type, subclass af.CrossSection, pass the
computed properties to super().__init__, and override stress.
Supports and loads
beam.fix(node)clamps a node (all 6 dofs);beam.pin(node)fixes its translations only.beam.add_load(loads)sets the nodal loads, shape(num_nodes, 6), ordered[Fx, Fy, Fz, Mx, My, Mz]in global coordinates. Distributed loads must be lumped to the nodes.beam.add_inertial_mass(masses)adds point masses at the nodes, shape(num_nodes,).af.Frame(..., acc=acc)applies a frame acceleration[ax, ay, az, αx, αy, αz], shape(6,), as inertial loadsM @ accon the structure and on the point masses. Use[0, 0, -9.81, 0, 0, 0]for gravity.
Joints and frames
af.Joint(members, nodes) rigidly connects one node of each member (negative indices
count from the end of the beam). af.Frame(beams, joints, acc) assembles everything.
A portal frame under gravity and a lateral load:
import numpy as np
import csdl_alpha as csdl
import aframe as af
recorder = csdl.Recorder(inline=True)
recorder.start()
n = 11
def tube_beam(name, start, end, z=False):
mesh = csdl.Variable(value=np.linspace(start, end, n))
cs = af.CSTube(radius=csdl.Variable(value=np.full(n - 1, 0.15)),
thickness=csdl.Variable(value=np.full(n - 1, 0.01)))
return af.Beam(name, mesh, E=69e9, G=26e9, density=2700, cs=cs, z=z)
left = tube_beam('left', [0, 0, 0], [0, 0, 5], z=True) # vertical members need z=True
right = tube_beam('right', [8, 0, 0], [8, 0, 5], z=True)
top = tube_beam('top', [0, 0, 5], [8, 0, 5])
left.fix(0)
right.fix(0)
loads = np.zeros((n, 6))
loads[n // 2, 1] = 5000 # 5 kN lateral load at midspan
top.add_load(csdl.Variable(value=loads))
joints = [af.Joint(members=[left, top], nodes=[-1, 0]), # rigid corners
af.Joint(members=[right, top], nodes=[-1, -1])]
gravity = csdl.Variable(value=np.array([0, 0, -9.81, 0, 0, 0]))
frame = af.Frame(beams=[left, right, top], joints=joints, acc=gravity)
frame.solve()
stress = frame.compute_stress()
cg, mass = frame.compute_mass_properties()
recorder.stop()
print('midspan displacement (m):', frame.displacement['top'].value[n // 2])
print('max stress (Pa):', max(s.value.max() for s in stress.values()))
print('mass (kg):', mass.value, 'cg (m):', cg.value)
A beam can belong to several frames, e.g. for different joint layouts.
Results
| Result | Description |
|---|---|
frame.displacement[name], frame.rotation[name] |
Nodal displacements and rotations of a beam, shape (num_nodes, 3) (after solve()) |
frame.compute_stress(beams=None) |
Dict of per-beam stresses (see the cross-section table); pass beams to skip beams without a stress model |
frame.compute_mass_properties() |
Center of gravity (3,) and total mass of the beams (point masses excluded) |
beam.mass, beam.cg |
Mass and center of gravity of a beam |
frame.U, frame.K, frame.M, frame.F |
Global displacement vector, stiffness and mass matrices, and load vector |
frame.node_dofs(beam) |
Global dof indices of a beam's nodes, shape (num_nodes, 6), e.g. for coupling to other solvers |
All results are CSDL variables: read their .value, differentiate them with
csdl.derivative, or use them as objectives and constraints.
Optimization
Set inputs as design variables and outputs as objectives or constraints, then hand the recorder to a CSDL simulator and an optimizer. A minimum-mass cantilever with stress and deflection limits, using modOpt and the JAX backend:
import numpy as np
import csdl_alpha as csdl
import aframe as af
from modopt import CSDLAlphaProblem, SLSQP
recorder = csdl.Recorder(inline=True)
recorder.start()
num_nodes = 21
mesh = np.zeros((num_nodes, 3))
mesh[:, 0] = np.linspace(0, 10, num_nodes)
radius = csdl.Variable(value=np.full(num_nodes - 1, 0.1), name='radius')
radius.set_as_design_variable(lower=0.02, upper=0.5, scaler=10)
thickness = csdl.Variable(value=np.full(num_nodes - 1, 0.005), name='thickness')
thickness.set_as_design_variable(lower=0.001, upper=0.02, scaler=100)
beam = af.Beam(name='cantilever', mesh=csdl.Variable(value=mesh), E=69e9, G=26e9,
density=2700, cs=af.CSTube(radius=radius, thickness=thickness))
beam.fix(0)
loads = np.zeros((num_nodes, 6))
loads[-1, 2] = -1000
beam.add_load(csdl.Variable(value=loads))
frame = af.Frame(beams=[beam])
frame.solve()
stress = frame.compute_stress()['cantilever']
stress.set_as_constraint(upper=100e6, scaler=1e-8) # 100 MPa allowable
tip_deflection = frame.displacement['cantilever'][-1, 2]
tip_deflection.set_as_constraint(lower=-0.2, scaler=10) # at most 0.2 m
beam.mass.set_as_objective(scaler=1e-2)
recorder.stop()
sim = csdl.experimental.JaxSimulator(recorder=recorder)
problem = CSDLAlphaProblem(problem_name='cantilever', simulator=sim)
optimizer = SLSQP(problem, solver_options={'maxiter': 200, 'ftol': 1e-8}, turn_off_outputs=True)
optimizer.solve()
optimizer.print_results()
csdl.experimental.PySimulator(recorder) runs the same model without JAX. See
examples/single_beam_opt.py for a complete script.
Plotting
With the plot extra, the PyVista helpers draw beams into a pyvista.Plotter. For
example, the deformed cantilever from the quick start, colored by stress:
import pyvista as pv
deformed = mesh + frame.displacement['cantilever'].value
plotter = pv.Plotter()
af.plot_mesh(plotter, deformed, cell_data=stress.value, line_width=10)
af.plot_points(plotter, deformed, point_size=15, color='black')
plotter.show()
af.plot_cyl and af.plot_box draw each element as a cylinder or box, and
aframe.utils.plot_matplotlib has Matplotlib equivalents.
Modeling assumptions
- Linear statics: small displacements and rotations, linear elastic material.
- Euler-Bernoulli elements (no shear deformation); cross-section properties are constant within an element.
- Joints are rigid and share all 6 dofs.
- Stresses are evaluated at element midpoints from the element internal loads.
- The global matrices are dense and solved with a dense direct solver, so time and memory grow quickly with model size; models of up to a few thousand dofs (several hundred nodes) solve in seconds.
Documentation
The API documentation is at https://lsdolab.github.io/aframe/. It is generated from the
docstrings with pdoc and republished on every push to main. To build it
locally:
pip install -e ".[plot,docs]"
python docs/make_docs.py
and open site/index.html.
Examples
The examples folder has complete scripts:
| Script | What it shows |
|---|---|
two_joined_beams.py |
Two tube beams connected by a joint |
bingran_cantilever_beam.py |
A tube cantilever with stresses and mass properties |
single_beam_opt.py |
Minimum-mass optimization with a displacement constraint (modOpt) |
NASA_LPC_beam.py, aurora_pav_vv.py |
Wing-box (CSBox) beams under spanwise load distributions |
1_element_beam_test.py |
Frames of single-element beams |
ucsd_lunar_lander/ |
A lunar lander frame with many joints, plotted with PyVista |
bunny/ |
A frame generated from the edges of an STL mesh |
License
This project is licensed under the terms of the MIT License.
1""" 2aframe: a differentiable linear 3D beam/frame solver written in CSDL. 3 4.. include:: ../README.md 5 :start-line: 1 6""" 7from aframe.core.cs import CrossSection, CSTube, CSCircle, CSEllipse, CSBox, CSBoxMarius 8from aframe.core.beam import Beam 9from aframe.core.joint import Joint 10from aframe.core.frame import Frame 11from aframe.utils.meshing import mesh_from_points_and_edges 12 13__version__ = '0.1.0' 14 15# the PyVista plotting helpers need the optional 'plot' dependencies and are 16# slow to import, so they are loaded on first use 17_PLOTTING = ('plot_box', 'plot_cyl', 'plot_mesh', 'plot_points') 18 19__all__ = ['CrossSection', 'CSTube', 'CSCircle', 'CSEllipse', 'CSBox', 'CSBoxMarius', 20 'Beam', 'Joint', 'Frame', 'mesh_from_points_and_edges', *_PLOTTING] 21 22 23def __getattr__(name): 24 if name in _PLOTTING: 25 try: 26 from aframe.utils import plot_pyvista 27 except ImportError as e: 28 raise ImportError(f"aframe.{name} needs pyvista and matplotlib: pip install aframe[plot]") from e 29 return getattr(plot_pyvista, name) 30 31 raise AttributeError(f"module 'aframe' has no attribute {name!r}") 32 33 34def __dir__(): 35 return sorted(set(globals()) | set(__all__))
class
CrossSection:
23class CrossSection: 24 """ 25 Base class of all cross-sections, also usable directly for a section with 26 known properties (it has no stress model). 27 28 Subclasses compute their properties and pass them to ``__init__``; those 29 with a stress model override ``stress``. 30 31 Parameters 32 ---------- 33 area, ix, iy, iz : csdl.Variable 34 Section properties, shape (num_elements,). 35 """ 36 37 def __init__(self, area:csdl.Variable, ix:csdl.Variable, iy:csdl.Variable, iz:csdl.Variable): 38 39 if not area.shape == ix.shape == iy.shape == iz.shape: 40 raise ValueError('section properties area, ix, iy and iz must have the same shape') 41 42 self.area = area 43 self.ix = ix 44 self.iy = iy 45 self.iz = iz 46 47 48 @property 49 def num_elements(self)->int: 50 return self.area.shape[0] 51 52 53 @staticmethod 54 def internal_loads(element_loads:csdl.Variable)->Tuple[csdl.Variable, ...]: 55 """ 56 internal loads at the element midspan from the local element end loads 57 [F_a, M_a, F_b, M_b], shape (num_elements, 12) 58 59 returns (N, V_y, V_z, T, M_y, M_z), each shape (num_elements,): axial 60 force (tension positive), shear forces, torque and bending moments. 61 The end loads act on the element, so the internal load is (end_b - end_a) / 2 62 (the forces and torque are constant along an element, the moments linear). 63 """ 64 return tuple((element_loads[:, 6 + i] - element_loads[:, i]) / 2 for i in range(6)) 65 66 67 def stress(self, element_loads:csdl.Variable)->csdl.Variable: 68 raise NotImplementedError(f'{type(self).__name__} has no stress model')
Base class of all cross-sections, also usable directly for a section with known properties (it has no stress model).
Subclasses compute their properties and pass them to __init__; those
with a stress model override stress.
Parameters
- area, ix, iy, iz (csdl.Variable): Section properties, shape (num_elements,).
CrossSection( area: csdl_alpha.src.graph.variable.Variable, ix: csdl_alpha.src.graph.variable.Variable, iy: csdl_alpha.src.graph.variable.Variable, iz: csdl_alpha.src.graph.variable.Variable)
37 def __init__(self, area:csdl.Variable, ix:csdl.Variable, iy:csdl.Variable, iz:csdl.Variable): 38 39 if not area.shape == ix.shape == iy.shape == iz.shape: 40 raise ValueError('section properties area, ix, iy and iz must have the same shape') 41 42 self.area = area 43 self.ix = ix 44 self.iy = iy 45 self.iz = iz
@staticmethod
def
internal_loads( element_loads: csdl_alpha.src.graph.variable.Variable) -> Tuple[csdl_alpha.src.graph.variable.Variable, ...]:
53 @staticmethod 54 def internal_loads(element_loads:csdl.Variable)->Tuple[csdl.Variable, ...]: 55 """ 56 internal loads at the element midspan from the local element end loads 57 [F_a, M_a, F_b, M_b], shape (num_elements, 12) 58 59 returns (N, V_y, V_z, T, M_y, M_z), each shape (num_elements,): axial 60 force (tension positive), shear forces, torque and bending moments. 61 The end loads act on the element, so the internal load is (end_b - end_a) / 2 62 (the forces and torque are constant along an element, the moments linear). 63 """ 64 return tuple((element_loads[:, 6 + i] - element_loads[:, i]) / 2 for i in range(6))
internal loads at the element midspan from the local element end loads [F_a, M_a, F_b, M_b], shape (num_elements, 12)
returns (N, V_y, V_z, T, M_y, M_z), each shape (num_elements,): axial force (tension positive), shear forces, torque and bending moments. The end loads act on the element, so the internal load is (end_b - end_a) / 2 (the forces and torque are constant along an element, the moments linear).
71class CSTube(CrossSection): 72 """ 73 Hollow circular tube. 74 75 Parameters 76 ---------- 77 radius : csdl.Variable 78 Outer radius, shape (num_elements,). 79 thickness : csdl.Variable 80 Wall thickness, shape (num_elements,). 81 """ 82 83 def __init__(self, radius:csdl.Variable, thickness:csdl.Variable): 84 85 if radius.shape != thickness.shape: 86 raise ValueError('tube radius and thickness must have the same shape') 87 88 self.radius = radius 89 self.thickness = thickness 90 self.inner_radius = radius - thickness 91 self.outer_radius = radius 92 self.precomp = np.pi * (self.outer_radius**4 - self.inner_radius**4) 93 94 super().__init__(area=np.pi * (self.outer_radius**2 - self.inner_radius**2), 95 ix=self.precomp / 2, 96 iy=self.precomp / 4, 97 iz=self.precomp / 4) 98 99 100 def stress(self, element_loads:csdl.Variable)->csdl.Variable: 101 """ 102 von Mises stress at the most stressed point of the outer fiber, 103 shape (num_elements,): |N| / A + M r / I combined with the torsional shear 104 """ 105 N, _, _, T, M_y, M_z = self.internal_loads(element_loads) 106 107 # eps keeps the absolute values/square roots differentiable at zero load 108 eps = 1E-12 109 axial_stress = (N**2 + eps) ** 0.5 / self.area 110 shear_stress = T * self.radius / self.ix 111 112 max_moment = (M_y**2 + M_z**2 + eps) ** 0.5 113 bending_stress = max_moment * self.radius / self.iy 114 115 normal_stress = axial_stress + bending_stress 116 117 return (normal_stress**2 + 3*shear_stress**2 + eps) ** 0.5
Hollow circular tube.
Parameters
- radius (csdl.Variable): Outer radius, shape (num_elements,).
- thickness (csdl.Variable): Wall thickness, shape (num_elements,).
CSTube( radius: csdl_alpha.src.graph.variable.Variable, thickness: csdl_alpha.src.graph.variable.Variable)
83 def __init__(self, radius:csdl.Variable, thickness:csdl.Variable): 84 85 if radius.shape != thickness.shape: 86 raise ValueError('tube radius and thickness must have the same shape') 87 88 self.radius = radius 89 self.thickness = thickness 90 self.inner_radius = radius - thickness 91 self.outer_radius = radius 92 self.precomp = np.pi * (self.outer_radius**4 - self.inner_radius**4) 93 94 super().__init__(area=np.pi * (self.outer_radius**2 - self.inner_radius**2), 95 ix=self.precomp / 2, 96 iy=self.precomp / 4, 97 iz=self.precomp / 4)
def
stress( self, element_loads: csdl_alpha.src.graph.variable.Variable) -> csdl_alpha.src.graph.variable.Variable:
100 def stress(self, element_loads:csdl.Variable)->csdl.Variable: 101 """ 102 von Mises stress at the most stressed point of the outer fiber, 103 shape (num_elements,): |N| / A + M r / I combined with the torsional shear 104 """ 105 N, _, _, T, M_y, M_z = self.internal_loads(element_loads) 106 107 # eps keeps the absolute values/square roots differentiable at zero load 108 eps = 1E-12 109 axial_stress = (N**2 + eps) ** 0.5 / self.area 110 shear_stress = T * self.radius / self.ix 111 112 max_moment = (M_y**2 + M_z**2 + eps) ** 0.5 113 bending_stress = max_moment * self.radius / self.iy 114 115 normal_stress = axial_stress + bending_stress 116 117 return (normal_stress**2 + 3*shear_stress**2 + eps) ** 0.5
von Mises stress at the most stressed point of the outer fiber, shape (num_elements,): |N| / A + M r / I combined with the torsional shear
120class CSCircle(CrossSection): 121 """ 122 Solid circular section (no stress model). 123 124 Parameters 125 ---------- 126 radius : csdl.Variable 127 Radius, shape (num_elements,). 128 """ 129 130 def __init__(self, radius:csdl.Variable): 131 132 self.radius = radius 133 self.precomp = np.pi * self.radius**4 134 135 super().__init__(area=np.pi * self.radius**2, 136 ix=(1 / 2) * self.precomp, 137 iy=(1 / 4) * self.precomp, 138 iz=(1 / 4) * self.precomp)
Solid circular section (no stress model).
Parameters
- radius (csdl.Variable): Radius, shape (num_elements,).
141class CSEllipse(CrossSection): 142 """ 143 Solid elliptical section (no stress model). 144 145 Parameters 146 ---------- 147 semi_major_axis : csdl.Variable 148 Semi-axis along the local z direction, shape (num_elements,). 149 semi_minor_axis : csdl.Variable 150 Semi-axis along the local y direction, shape (num_elements,). 151 """ 152 153 def __init__(self, semi_major_axis:csdl.Variable, semi_minor_axis:csdl.Variable): 154 155 a = self.semi_major_axis = semi_major_axis 156 b = self.semi_minor_axis = semi_minor_axis 157 158 area = np.pi * a * b 159 # approximate torsion constant 160 beta = 1 / ((1 + (b / a)**2)**0.5) 161 ix = (np.pi / 2) * a * b**3 * beta 162 iy = np.pi / 4 * a * b**3 163 iz = np.pi / 4 * a**3 * b 164 165 super().__init__(area=area, ix=ix, iy=iy, iz=iz)
Solid elliptical section (no stress model).
Parameters
- semi_major_axis (csdl.Variable): Semi-axis along the local z direction, shape (num_elements,).
- semi_minor_axis (csdl.Variable): Semi-axis along the local y direction, shape (num_elements,).
CSEllipse( semi_major_axis: csdl_alpha.src.graph.variable.Variable, semi_minor_axis: csdl_alpha.src.graph.variable.Variable)
153 def __init__(self, semi_major_axis:csdl.Variable, semi_minor_axis:csdl.Variable): 154 155 a = self.semi_major_axis = semi_major_axis 156 b = self.semi_minor_axis = semi_minor_axis 157 158 area = np.pi * a * b 159 # approximate torsion constant 160 beta = 1 / ((1 + (b / a)**2)**0.5) 161 ix = (np.pi / 2) * a * b**3 * beta 162 iy = np.pi / 4 * a * b**3 163 iz = np.pi / 4 * a**3 * b 164 165 super().__init__(area=area, ix=ix, iy=iy, iz=iz)
168class CSBox(CrossSection): 169 """ 170 Thin-walled rectangular box (e.g. a wing box), with separate top/bottom 171 skin and web thicknesses. Use ``CSBox.from_properties`` to supply the 172 section properties instead of the thin-walled formulas. 173 174 Parameters 175 ---------- 176 ttop, tbot : csdl.Variable 177 Top and bottom skin thicknesses, shape (num_elements,). 178 tweb : csdl.Variable 179 Web (side wall) thickness, shape (num_elements,). 180 height, width : csdl.Variable 181 Outer height (local z) and width (local y), shape (num_elements,). 182 """ 183 184 def __init__(self, 185 ttop:csdl.Variable, 186 tbot:csdl.Variable, 187 tweb:csdl.Variable, 188 height:csdl.Variable, 189 width:csdl.Variable): 190 191 self._set_dimensions(ttop, tbot, tweb, height, width) 192 193 self.neutral_axis = self._neutral_axis() 194 area = self._area() 195 iy = self._iy() 196 iz = self._iz() 197 # thin-walled approximation 198 super().__init__(area=area, ix=iy + iz, iy=iy, iz=iz) 199 200 201 @classmethod 202 def from_properties(cls, 203 ttop:csdl.Variable, 204 tbot:csdl.Variable, 205 tweb:csdl.Variable, 206 height:csdl.Variable, 207 width:csdl.Variable, 208 area:csdl.Variable, 209 ix:csdl.Variable, 210 iy:csdl.Variable, 211 iz:csdl.Variable)->'CSBox': 212 """a box with user-supplied section properties; the dimensions are used for stress recovery""" 213 box = cls.__new__(cls) 214 box._set_dimensions(ttop, tbot, tweb, height, width) 215 CrossSection.__init__(box, area, ix, iy, iz) 216 return box 217 218 219 def _set_dimensions(self, ttop, tbot, tweb, height, width): 220 if not ttop.shape == tbot.shape == tweb.shape == height.shape == width.shape: 221 raise ValueError('box beam ttop, tbot, tweb, height and width must have the same shape') 222 223 self.ttop = ttop 224 self.tbot = tbot 225 self.tweb = tweb 226 self.height = height 227 self.width = width 228 229 230 def _neutral_axis(self)->tuple: 231 """(y, z) offset of the neutral axis from the box center""" 232 Atop = self.ttop * (self.width - 2 * self.tweb) 233 Abot = self.tbot * (self.width - 2 * self.tweb) 234 235 centroid_z = (Atop - Abot) / self.width 236 centroid_y = 0 237 238 return (centroid_y, centroid_z) 239 240 241 def _area(self)->csdl.Variable: 242 w_i = self.width - 2 * self.tweb 243 h_i = self.height - self.ttop - self.tbot 244 return self.width * self.height - w_i * h_i 245 246 247 def _iy(self)->csdl.Variable: 248 """parallel-axis sum of the skins and webs about the neutral axis""" 249 cz_neutral = self.neutral_axis[1] 250 251 # top skin 252 area_top = self.ttop * self.width 253 iy_top = self.ttop ** 3 * self.width / 12 254 d_top = self.height / 2 - self.ttop / 2 - cz_neutral 255 iy_top_centroid = iy_top + area_top * d_top**2 256 257 # bottom skin 258 area_bot = self.tbot * self.width 259 iy_bot = self.tbot ** 3 * self.width / 12 260 d_bot = -self.height / 2 + self.tbot / 2 - cz_neutral 261 iy_bot_centroid = iy_bot + area_bot * d_bot**2 262 263 # front/rear webs (centered at z = 0) 264 area_web = self.tweb * (self.height - self.ttop - self.tbot) 265 iy_web = (self.height - self.ttop - self.tbot) ** 3 * self.tweb / 12 266 iy_rear_centroid = iy_web + area_web * cz_neutral**2 267 iy_front_centroid = iy_web + area_web * cz_neutral**2 268 269 return iy_top_centroid + iy_bot_centroid + iy_rear_centroid + iy_front_centroid 270 271 272 def _iz(self)->csdl.Variable: 273 """parallel-axis sum of the skins and webs about the neutral axis""" 274 cy_neutral = self.neutral_axis[0] 275 276 # top/bottom skins (centered at y = 0) 277 area_top = self.ttop * self.width 278 iz_top = self.ttop * self.width**3 / 12 279 area_bot = self.tbot * self.width 280 iz_bot = self.tbot * self.width**3 / 12 281 282 iz_top_centroid = iz_top + area_top * cy_neutral**2 283 iz_bot_centroid = iz_bot + area_bot * cy_neutral**2 284 285 # front/rear webs 286 area_web = self.tweb * (self.height - self.ttop - self.tbot) 287 iz_web = (self.height - self.ttop - self.tbot) * self.tweb**3 / 12 288 289 d_front = -self.width / 2 + self.tweb / 2 - cy_neutral 290 d_rear = self.width / 2 - self.tweb / 2 - cy_neutral 291 292 iz_front_centroid = iz_web + area_web * d_front**2 293 iz_rear_centroid = iz_web + area_web * d_rear**2 294 295 return iz_top_centroid + iz_bot_centroid + iz_front_centroid + iz_rear_centroid 296 297 298 def stress(self, element_loads:csdl.Variable)->csdl.Variable: 299 """ 300 von Mises stress at 5 points of the section, shape (num_elements, 5) 301 302 0-----------------1 303 | | 304 4 | 305 | | 306 3-----------------2 307 308 points 0-3 are the corners (axial + bending + torsion); point 4 is the 309 web mid-height, which adds the vertical shear from V_z 310 """ 311 N, _, V_z, T, M_y, M_z = self.internal_loads(element_loads) 312 313 # the axial stress is common to all stress evaluation points 314 axial_stress = N / self.area 315 316 def von_mises(z, y, shear=None): 317 p = (z**2 + y**2)**0.5 318 torsional_stress = T * p / self.ix 319 if shear is not None: 320 torsional_stress = torsional_stress + shear 321 normal_stress = axial_stress + M_y * y / self.iy + M_z * z / self.iz 322 # 1E-8 keeps the square root differentiable at zero stress 323 return (normal_stress**2 + 3*torsional_stress**2 + 1E-8)**0.5 324 325 corners = [(-self.width / 2, self.height / 2), 326 (self.width / 2, self.height / 2), 327 (self.width / 2, -self.height / 2), 328 (-self.width / 2, -self.height / 2)] 329 points = [von_mises(z, y) for z, y in corners] 330 331 # web mid-height: add the vertical shear V_z Q / (I t), with an 332 # approximate first moment of area Q and both webs (t = 2 tweb) 333 tcap = (self.ttop + self.tbot) / 2 334 Q = self.width * tcap * (self.height / 2) + 2 * (self.height / 2) * self.tweb * (self.height / 4) 335 shear_stress = V_z * Q / (self.iy * 2 * self.tweb) 336 points.append(von_mises(-self.width / 2, 0, shear_stress)) 337 338 return csdl.transpose(csdl.vstack(points))
Thin-walled rectangular box (e.g. a wing box), with separate top/bottom
skin and web thicknesses. Use CSBox.from_properties to supply the
section properties instead of the thin-walled formulas.
Parameters
- ttop, tbot (csdl.Variable): Top and bottom skin thicknesses, shape (num_elements,).
- tweb (csdl.Variable): Web (side wall) thickness, shape (num_elements,).
- height, width (csdl.Variable): Outer height (local z) and width (local y), shape (num_elements,).
CSBox( ttop: csdl_alpha.src.graph.variable.Variable, tbot: csdl_alpha.src.graph.variable.Variable, tweb: csdl_alpha.src.graph.variable.Variable, height: csdl_alpha.src.graph.variable.Variable, width: csdl_alpha.src.graph.variable.Variable)
184 def __init__(self, 185 ttop:csdl.Variable, 186 tbot:csdl.Variable, 187 tweb:csdl.Variable, 188 height:csdl.Variable, 189 width:csdl.Variable): 190 191 self._set_dimensions(ttop, tbot, tweb, height, width) 192 193 self.neutral_axis = self._neutral_axis() 194 area = self._area() 195 iy = self._iy() 196 iz = self._iz() 197 # thin-walled approximation 198 super().__init__(area=area, ix=iy + iz, iy=iy, iz=iz)
@classmethod
def
from_properties( cls, ttop: csdl_alpha.src.graph.variable.Variable, tbot: csdl_alpha.src.graph.variable.Variable, tweb: csdl_alpha.src.graph.variable.Variable, height: csdl_alpha.src.graph.variable.Variable, width: csdl_alpha.src.graph.variable.Variable, area: csdl_alpha.src.graph.variable.Variable, ix: csdl_alpha.src.graph.variable.Variable, iy: csdl_alpha.src.graph.variable.Variable, iz: csdl_alpha.src.graph.variable.Variable) -> CSBox:
201 @classmethod 202 def from_properties(cls, 203 ttop:csdl.Variable, 204 tbot:csdl.Variable, 205 tweb:csdl.Variable, 206 height:csdl.Variable, 207 width:csdl.Variable, 208 area:csdl.Variable, 209 ix:csdl.Variable, 210 iy:csdl.Variable, 211 iz:csdl.Variable)->'CSBox': 212 """a box with user-supplied section properties; the dimensions are used for stress recovery""" 213 box = cls.__new__(cls) 214 box._set_dimensions(ttop, tbot, tweb, height, width) 215 CrossSection.__init__(box, area, ix, iy, iz) 216 return box
a box with user-supplied section properties; the dimensions are used for stress recovery
def
stress( self, element_loads: csdl_alpha.src.graph.variable.Variable) -> csdl_alpha.src.graph.variable.Variable:
298 def stress(self, element_loads:csdl.Variable)->csdl.Variable: 299 """ 300 von Mises stress at 5 points of the section, shape (num_elements, 5) 301 302 0-----------------1 303 | | 304 4 | 305 | | 306 3-----------------2 307 308 points 0-3 are the corners (axial + bending + torsion); point 4 is the 309 web mid-height, which adds the vertical shear from V_z 310 """ 311 N, _, V_z, T, M_y, M_z = self.internal_loads(element_loads) 312 313 # the axial stress is common to all stress evaluation points 314 axial_stress = N / self.area 315 316 def von_mises(z, y, shear=None): 317 p = (z**2 + y**2)**0.5 318 torsional_stress = T * p / self.ix 319 if shear is not None: 320 torsional_stress = torsional_stress + shear 321 normal_stress = axial_stress + M_y * y / self.iy + M_z * z / self.iz 322 # 1E-8 keeps the square root differentiable at zero stress 323 return (normal_stress**2 + 3*torsional_stress**2 + 1E-8)**0.5 324 325 corners = [(-self.width / 2, self.height / 2), 326 (self.width / 2, self.height / 2), 327 (self.width / 2, -self.height / 2), 328 (-self.width / 2, -self.height / 2)] 329 points = [von_mises(z, y) for z, y in corners] 330 331 # web mid-height: add the vertical shear V_z Q / (I t), with an 332 # approximate first moment of area Q and both webs (t = 2 tweb) 333 tcap = (self.ttop + self.tbot) / 2 334 Q = self.width * tcap * (self.height / 2) + 2 * (self.height / 2) * self.tweb * (self.height / 4) 335 shear_stress = V_z * Q / (self.iy * 2 * self.tweb) 336 points.append(von_mises(-self.width / 2, 0, shear_stress)) 337 338 return csdl.transpose(csdl.vstack(points))
von Mises stress at 5 points of the section, shape (num_elements, 5)
0-----------------1
| |
4 |
| |
3-----------------2
points 0-3 are the corners (axial + bending + torsion); point 4 is the web mid-height, which adds the vertical shear from V_z
341class CSBoxMarius(CSBox): 342 """ 343 Box with user-supplied section properties. 344 345 Deprecated: use ``CSBox.from_properties``, which takes the same arguments. 346 """ 347 348 def __init__(self, 349 ttop:csdl.Variable, 350 tbot:csdl.Variable, 351 tweb:csdl.Variable, 352 height:csdl.Variable, 353 width:csdl.Variable, 354 area:csdl.Variable, 355 ix:csdl.Variable, 356 iy:csdl.Variable, 357 iz:csdl.Variable): 358 359 self._set_dimensions(ttop, tbot, tweb, height, width) 360 CrossSection.__init__(self, area, ix, iy, iz)
Box with user-supplied section properties.
Deprecated: use CSBox.from_properties, which takes the same arguments.
CSBoxMarius( ttop: csdl_alpha.src.graph.variable.Variable, tbot: csdl_alpha.src.graph.variable.Variable, tweb: csdl_alpha.src.graph.variable.Variable, height: csdl_alpha.src.graph.variable.Variable, width: csdl_alpha.src.graph.variable.Variable, area: csdl_alpha.src.graph.variable.Variable, ix: csdl_alpha.src.graph.variable.Variable, iy: csdl_alpha.src.graph.variable.Variable, iz: csdl_alpha.src.graph.variable.Variable)
348 def __init__(self, 349 ttop:csdl.Variable, 350 tbot:csdl.Variable, 351 tweb:csdl.Variable, 352 height:csdl.Variable, 353 width:csdl.Variable, 354 area:csdl.Variable, 355 ix:csdl.Variable, 356 iy:csdl.Variable, 357 iz:csdl.Variable): 358 359 self._set_dimensions(ttop, tbot, tweb, height, width) 360 CrossSection.__init__(self, area, ix, iy, iz)
class
Beam:
14class Beam: 15 """ 16 A beam discretized into num_nodes - 1 two-node Euler-Bernoulli elements 17 with 6 dofs per node (3 translations, 3 rotations). 18 19 Element matrices are built in local coordinates (x along the element) and 20 rotated to global coordinates. The mesh, material and cross-section are 21 csdl Variables/constants, so everything is differentiable. 22 23 Parameters 24 ---------- 25 name : str 26 Unique name, used as the key of the per-beam results in Frame. 27 mesh : csdl.Variable 28 Node coordinates, shape (num_nodes, 3). 29 E, G : float 30 Young's modulus and shear modulus. 31 density : float 32 Material density. 33 cs : CrossSection 34 Cross-section (CSTube, CSBox, ...) with per-element properties, 35 shape (num_nodes - 1,). 36 z : bool, optional 37 Set True for a beam along the global z axis, where the default 38 local-axis construction is singular. 39 """ 40 41 def __init__(self, 42 name:str, 43 mesh:csdl.Variable, 44 E:float, 45 G:float, 46 density:float, 47 cs:CrossSection, 48 z:bool = False): 49 50 self.name = name 51 self.mesh = mesh 52 self.E = E 53 self.G = G 54 self.density = density 55 self.cs = cs 56 self.z = z 57 58 self.num_nodes = mesh.shape[0] 59 self.num_elements = self.num_nodes - 1 60 61 if cs.area.shape != (self.num_elements,): 62 raise ValueError(f'the cross-section of beam {name!r} has properties of shape {cs.area.shape}, ' 63 f'expected (num_nodes - 1,) = ({self.num_elements},)') 64 65 self.loads = None 66 self.extra_inertial_mass = None 67 self.fixed_boundary_conditions: List[int] = [] 68 self.pinned_boundary_conditions: List[int] = [] 69 70 # element geometry and direction cosines 71 self.lengths, self.ll, self.mm, self.nn, self.D = self._lengths(mesh) 72 73 # element matrices in local and global coordinates 74 self.local_stiffness = self._local_stiffness_matrices() 75 self.local_mass = self._local_mass_matrices() 76 self.transforms = self._vectorized_transforms() 77 self.transformed_stiffness = self._transform_stiffness_matrices() 78 self.transformed_mass = self._transform_mass_matrices() 79 80 # mass properties 81 element_masses = self.cs.area * self.lengths * self.density 82 self.mass = csdl.sum(element_masses) 83 84 cg2 = (self.mesh[1:, :] + self.mesh[:-1, :]) / 2 85 element_masses_expanded = csdl.expand(element_masses, (self.num_elements, 3), action='i->ij') 86 self.rmvec = csdl.sum(cg2 * element_masses_expanded, axes=(0,)) 87 self.cg = self.rmvec / self.mass 88 89 90 def fix(self, node:int)->None: 91 """clamp a node (all 6 dofs)""" 92 if node < 0 or node > self.num_nodes - 1: 93 raise ValueError('fixed nodes must be between 0 and num_nodes - 1') 94 95 if node not in self.fixed_boundary_conditions and node not in self.pinned_boundary_conditions: 96 self.fixed_boundary_conditions.append(node) 97 98 99 def pin(self, node:int)->None: 100 """pin a node (the 3 translations; rotations stay free)""" 101 if node < 0 or node > self.num_nodes - 1: 102 raise ValueError('pinned nodes must be between 0 and num_nodes - 1') 103 104 if node not in self.pinned_boundary_conditions and node not in self.fixed_boundary_conditions: 105 self.pinned_boundary_conditions.append(node) 106 107 108 def add_inertial_mass(self, mass:csdl.Variable)->None: 109 """ 110 add point masses at the nodes, shape (num_nodes,); they are resolved 111 as inertial loads when the Frame has an acceleration 112 """ 113 if mass.shape != (self.num_nodes,): 114 raise ValueError('inertial mass must have shape (num_beam_nodes,)') 115 116 self.extra_inertial_mass = mass 117 118 119 def add_load(self, load:csdl.Variable)->None: 120 """set the nodal loads [Fx, Fy, Fz, Mx, My, Mz], shape (num_nodes, 6)""" 121 if load.shape != (self.num_nodes, 6): 122 raise ValueError('load must have shape (num_beam_nodes, 6)') 123 124 self.loads = load 125 126 127 def _lengths(self, mesh:csdl.Variable)->tuple: 128 """ 129 element lengths and direction cosines (ll, mm, nn) of the element axes, 130 plus D = sqrt(ll^2 + mm^2) used by the transforms 131 """ 132 diffs = mesh[1:] - mesh[:-1] 133 lengths = csdl.norm(diffs, axes=(1,)) 134 exl = csdl.expand(lengths, (self.num_elements, 3), action='i->ij') 135 cp = diffs / exl 136 137 ll = cp[:, 0] 138 mm = cp[:, 1] 139 nn = cp[:, 2] 140 D = (ll**2 + mm**2)**0.5 141 142 return lengths, ll, mm, nn, D 143 144 145 def _local_stiffness_matrices(self)->csdl.Variable: 146 """Euler-Bernoulli element stiffness matrices, shape (num_elements, 12, 12)""" 147 A = self.cs.area 148 E, G = self.E, self.G 149 Iz = self.cs.iz 150 Iy = self.cs.iy 151 J = self.cs.ix 152 L = self.lengths 153 154 AEL = A*E/L 155 nAEL = -AEL 156 GJL = G*J/L 157 nGJL = -GJL 158 159 # bending about local z (v, theta_z) 160 EIzL = E*Iz/L 161 EIzL2 = EIzL/L 162 EIzL3 = EIzL2/L 163 EIzL312 = 12*EIzL3 164 nEIzL312 = -EIzL312 165 EIzL26 = 6*EIzL2 166 nEIzL26 = -EIzL26 167 EIzL4 = 4*EIzL 168 EIzL_2 = 2*EIzL 169 170 # bending about local y (w, theta_y) 171 EIyL = E*Iy/L 172 EIyL2 = EIyL/L 173 EIyL3 = EIyL2/L 174 EIyL26 = 6*EIyL2 175 nEIyL26 = -EIyL26 176 EIyL312 = 12*EIyL3 177 nEIyL312 = -EIyL312 178 EIyL4 = 4*EIyL 179 EIyL_2 = 2*EIyL 180 181 # dof order per node: u, v, w, theta_x, theta_y, theta_z 182 diag = [(0, AEL), (1, EIzL312), (2, EIyL312), (3, GJL), (4, EIyL4), (5, EIzL4), 183 (6, AEL), (7, EIzL312), (8, EIyL312), (9, GJL), (10, EIyL4), (11, EIzL4)] 184 185 # upper triangle; the lower triangle is mirrored 186 off_diag = [(1, 5, EIzL26), (2, 4, nEIyL26), (0, 6, nAEL), (1, 7, nEIzL312), 187 (1, 11, EIzL26), (2, 8, nEIyL312), (2, 10, nEIyL26), (3, 9, nGJL), 188 (4, 8, EIyL26), (4, 10, EIyL_2), (5, 7, nEIzL26), (5, 11, EIzL_2), 189 (7, 11, nEIzL26), (8, 10, EIyL26)] 190 191 return self._symmetric_matrices(diag, off_diag) 192 193 194 def _local_mass_matrices(self)->csdl.Variable: 195 """consistent element mass matrices, shape (num_elements, 12, 12)""" 196 A = self.cs.area 197 rho = self.density 198 J = self.cs.ix 199 L = self.lengths 200 201 aa = L / 2 202 aa2 = aa**2 203 coef = rho * A * aa / 105 204 coef70 = coef * 70 205 coef78 = coef * 78 206 coef35 = coef * 35 207 ncoef35 = -coef35 208 coef27 = coef * 27 209 coef22aa = coef * 22 * aa 210 ncoef22aa = -coef22aa 211 coef13aa = coef * 13 * aa 212 ncoef13aa = -coef13aa 213 coef8aa2 = coef * 8 * aa2 214 ncoef6aa2 = -coef * 6 * aa2 215 # torsional inertia uses the polar radius of gyration squared 216 rx2 = J / A 217 coef70rx2 = coef70 * rx2 218 ncoef35rx2 = ncoef35 * rx2 219 220 diag = [(0, coef70), ([1, 2, 7, 8], coef78), (3, coef70rx2), 221 ([4, 5, 10, 11], coef8aa2), (6, coef70), (9, coef70rx2)] 222 223 # upper triangle; the lower triangle is mirrored 224 off_diag = [(2, 4, ncoef22aa), (1, 5, coef22aa), (0, 6, coef35), (1, 7, coef27), 225 (5, 7, coef13aa), (2, 8, coef27), (4, 8, ncoef13aa), (3, 9, ncoef35rx2), 226 (2, 10, coef13aa), (4, 10, ncoef6aa2), (8, 10, coef22aa), (1, 11, ncoef13aa), 227 (5, 11, ncoef6aa2), (7, 11, ncoef22aa)] 228 229 return self._symmetric_matrices(diag, off_diag) 230 231 232 def _symmetric_matrices(self, diag:list, off_diag:list)->csdl.Variable: 233 """ 234 build symmetric (num_elements, 12, 12) matrices from per-element values 235 diag: (index or list of indices, value); off_diag: (row, col, value), upper triangle 236 """ 237 n = self.num_elements 238 239 diag_matrices = csdl.Variable(value=np.zeros((n, 12, 12))) 240 for index, value in diag: 241 if isinstance(index, list): 242 value = value.expand((n, len(index)), action='i->ij') 243 diag_matrices = diag_matrices.set(csdl.slice[:, index, index], value) 244 245 upper = csdl.Variable(value=np.zeros((n, 12, 12))) 246 for row, col, value in off_diag: 247 upper = upper.set(csdl.slice[:, row, col], value) 248 249 return diag_matrices + upper + csdl.einsum(upper, action='ijk->ikj') 250 251 252 def _vectorized_transforms(self)->csdl.Variable: 253 """ 254 global-to-local rotation matrices T, shape (num_elements, 12, 12): 255 the same 3x3 direction-cosine block on each of the 4 diagonal blocks 256 """ 257 n = self.num_elements 258 259 if self.z: 260 # element along global z: local x = global z 261 zeros = csdl.Variable(value=np.zeros((n,))) 262 ones = csdl.Variable(value=np.ones((n,))) 263 block = [[zeros, zeros, ones], 264 [zeros, ones, None], 265 [-ones, zeros, zeros]] 266 else: 267 ll, mm, nn, D = self.ll, self.mm, self.nn, self.D 268 nmmD = -mm / D 269 llD = ll / D 270 block = [[ll, mm, nn], 271 [nmmD, llD, None], 272 [-nn * llD, nn * nmmD, D]] 273 274 T = csdl.Variable(value=np.zeros((n, 12, 12))) 275 for r in range(3): 276 for c in range(3): 277 if block[r][c] is not None: 278 # entry (r, c) of all 4 diagonal blocks 279 rows = [r, r + 3, r + 6, r + 9] 280 cols = [c, c + 3, c + 6, c + 9] 281 T = T.set(csdl.slice[:, rows, cols], block[r][c].expand((n, 4), action='i->ij')) 282 283 # kept for backward compatibility 284 self.transformations_bookshelf = T 285 286 return T 287 288 289 def _transform(self, local_matrices:csdl.Variable)->csdl.Variable: 290 """rotate local element matrices to global coordinates: T^T K T""" 291 transforms = self.transforms 292 T_transpose = csdl.einsum(transforms, action='ijk->ikj') 293 T_transpose_K = csdl.einsum(T_transpose, local_matrices, action='ijk,ikl->ijl') 294 return csdl.einsum(T_transpose_K, transforms, action='ijk,ikl->ijl') 295 296 297 def _transform_stiffness_matrices(self)->csdl.Variable: 298 """element stiffness matrices in global coordinates""" 299 return self._transform(self.local_stiffness) 300 301 302 def _transform_mass_matrices(self)->csdl.Variable: 303 """element mass matrices in global coordinates""" 304 return self._transform(self.local_mass) 305 306 307 def element_loads(self, element_displacements:csdl.Variable)->csdl.Variable: 308 """ 309 local element end loads K_local T u_e, shape (num_elements, 12), from the 310 global element displacements u_e = [u_a, u_b], shape (num_elements, 12) 311 """ 312 local_displacements = csdl.einsum(self.transforms, element_displacements, action='ijk,ik->ij') 313 314 return csdl.einsum(self.local_stiffness, local_displacements, action='ijk,ik->ij')
A beam discretized into num_nodes - 1 two-node Euler-Bernoulli elements with 6 dofs per node (3 translations, 3 rotations).
Element matrices are built in local coordinates (x along the element) and rotated to global coordinates. The mesh, material and cross-section are csdl Variables/constants, so everything is differentiable.
Parameters
- name (str): Unique name, used as the key of the per-beam results in Frame.
- mesh (csdl.Variable): Node coordinates, shape (num_nodes, 3).
- E, G (float): Young's modulus and shear modulus.
- density (float): Material density.
- cs (CrossSection): Cross-section (CSTube, CSBox, ...) with per-element properties, shape (num_nodes - 1,).
- z (bool, optional): Set True for a beam along the global z axis, where the default local-axis construction is singular.
Beam( name: str, mesh: csdl_alpha.src.graph.variable.Variable, E: float, G: float, density: float, cs: CrossSection, z: bool = False)
41 def __init__(self, 42 name:str, 43 mesh:csdl.Variable, 44 E:float, 45 G:float, 46 density:float, 47 cs:CrossSection, 48 z:bool = False): 49 50 self.name = name 51 self.mesh = mesh 52 self.E = E 53 self.G = G 54 self.density = density 55 self.cs = cs 56 self.z = z 57 58 self.num_nodes = mesh.shape[0] 59 self.num_elements = self.num_nodes - 1 60 61 if cs.area.shape != (self.num_elements,): 62 raise ValueError(f'the cross-section of beam {name!r} has properties of shape {cs.area.shape}, ' 63 f'expected (num_nodes - 1,) = ({self.num_elements},)') 64 65 self.loads = None 66 self.extra_inertial_mass = None 67 self.fixed_boundary_conditions: List[int] = [] 68 self.pinned_boundary_conditions: List[int] = [] 69 70 # element geometry and direction cosines 71 self.lengths, self.ll, self.mm, self.nn, self.D = self._lengths(mesh) 72 73 # element matrices in local and global coordinates 74 self.local_stiffness = self._local_stiffness_matrices() 75 self.local_mass = self._local_mass_matrices() 76 self.transforms = self._vectorized_transforms() 77 self.transformed_stiffness = self._transform_stiffness_matrices() 78 self.transformed_mass = self._transform_mass_matrices() 79 80 # mass properties 81 element_masses = self.cs.area * self.lengths * self.density 82 self.mass = csdl.sum(element_masses) 83 84 cg2 = (self.mesh[1:, :] + self.mesh[:-1, :]) / 2 85 element_masses_expanded = csdl.expand(element_masses, (self.num_elements, 3), action='i->ij') 86 self.rmvec = csdl.sum(cg2 * element_masses_expanded, axes=(0,)) 87 self.cg = self.rmvec / self.mass
def
fix(self, node: int) -> None:
90 def fix(self, node:int)->None: 91 """clamp a node (all 6 dofs)""" 92 if node < 0 or node > self.num_nodes - 1: 93 raise ValueError('fixed nodes must be between 0 and num_nodes - 1') 94 95 if node not in self.fixed_boundary_conditions and node not in self.pinned_boundary_conditions: 96 self.fixed_boundary_conditions.append(node)
clamp a node (all 6 dofs)
def
pin(self, node: int) -> None:
99 def pin(self, node:int)->None: 100 """pin a node (the 3 translations; rotations stay free)""" 101 if node < 0 or node > self.num_nodes - 1: 102 raise ValueError('pinned nodes must be between 0 and num_nodes - 1') 103 104 if node not in self.pinned_boundary_conditions and node not in self.fixed_boundary_conditions: 105 self.pinned_boundary_conditions.append(node)
pin a node (the 3 translations; rotations stay free)
def
add_inertial_mass(self, mass: csdl_alpha.src.graph.variable.Variable) -> None:
108 def add_inertial_mass(self, mass:csdl.Variable)->None: 109 """ 110 add point masses at the nodes, shape (num_nodes,); they are resolved 111 as inertial loads when the Frame has an acceleration 112 """ 113 if mass.shape != (self.num_nodes,): 114 raise ValueError('inertial mass must have shape (num_beam_nodes,)') 115 116 self.extra_inertial_mass = mass
add point masses at the nodes, shape (num_nodes,); they are resolved as inertial loads when the Frame has an acceleration
def
add_load(self, load: csdl_alpha.src.graph.variable.Variable) -> None:
119 def add_load(self, load:csdl.Variable)->None: 120 """set the nodal loads [Fx, Fy, Fz, Mx, My, Mz], shape (num_nodes, 6)""" 121 if load.shape != (self.num_nodes, 6): 122 raise ValueError('load must have shape (num_beam_nodes, 6)') 123 124 self.loads = load
set the nodal loads [Fx, Fy, Fz, Mx, My, Mz], shape (num_nodes, 6)
def
element_loads( self, element_displacements: csdl_alpha.src.graph.variable.Variable) -> csdl_alpha.src.graph.variable.Variable:
307 def element_loads(self, element_displacements:csdl.Variable)->csdl.Variable: 308 """ 309 local element end loads K_local T u_e, shape (num_elements, 12), from the 310 global element displacements u_e = [u_a, u_b], shape (num_elements, 12) 311 """ 312 local_displacements = csdl.einsum(self.transforms, element_displacements, action='ijk,ik->ij') 313 314 return csdl.einsum(self.local_stiffness, local_displacements, action='ijk,ik->ij')
local element end loads K_local T u_e, shape (num_elements, 12), from the global element displacements u_e = [u_a, u_b], shape (num_elements, 12)
class
Joint:
12class Joint: 13 """ 14 Rigidly connect beams by merging one node of each into a single frame node 15 (all 6 dofs are shared). 16 17 Parameters 18 ---------- 19 members : list of Beam 20 The beams to connect. 21 nodes : list of int 22 The node index in each member that is connected, same length as members 23 (negative indices count from the end of the beam). 24 """ 25 26 def __init__(self, members:List[Beam], nodes:List[int])->None: 27 28 if len(members) != len(nodes): 29 raise ValueError('a joint needs one node index per member') 30 31 for member, node in zip(members, nodes): 32 if not -member.num_nodes <= node < member.num_nodes: 33 raise ValueError(f'joint node {node} is out of range for beam {member.name!r} ' 34 f'with {member.num_nodes} nodes') 35 36 self.members = members 37 self.nodes = [int(node) % member.num_nodes for member, node in zip(members, nodes)]
Rigidly connect beams by merging one node of each into a single frame node (all 6 dofs are shared).
Parameters
- members (list of Beam): The beams to connect.
- nodes (list of int): The node index in each member that is connected, same length as members (negative indices count from the end of the beam).
Joint(members: List[Beam], nodes: List[int])
26 def __init__(self, members:List[Beam], nodes:List[int])->None: 27 28 if len(members) != len(nodes): 29 raise ValueError('a joint needs one node index per member') 30 31 for member, node in zip(members, nodes): 32 if not -member.num_nodes <= node < member.num_nodes: 33 raise ValueError(f'joint node {node} is out of range for beam {member.name!r} ' 34 f'with {member.num_nodes} nodes') 35 36 self.members = members 37 self.nodes = [int(node) % member.num_nodes for member, node in zip(members, nodes)]
class
Frame:
48class Frame: 49 """ 50 A linear static frame model built from beams and joints. 51 52 Building the Frame numbers the frame nodes and assembles the global 53 stiffness and mass matrices and the load vector (with boundary conditions 54 applied); call ``solve()`` for the displacements and ``compute_stress()`` 55 for the stresses. The dof numbering is owned by the Frame, so a Beam can 56 be used in several Frames. 57 58 Parameters 59 ---------- 60 beams : list of Beam 61 The beams (unique objects with unique names), with loads and boundary 62 conditions already added. 63 joints : list of Joint, optional 64 Rigid connections between beams of this frame. 65 acc : csdl.Variable, optional 66 Frame acceleration [ax, ay, az, alpha_x, alpha_y, alpha_z], shape (6,), 67 applied as inertial loads M @ acc on the structure and extra masses. 68 69 Attributes 70 ---------- 71 K, M : csdl.Variable 72 Global stiffness and mass matrices, shape (dim, dim), dim = 6 * num. 73 F : csdl.Variable 74 Global load vector, shape (dim,). 75 U : csdl.Variable 76 Global displacement vector, set by ``solve()``. 77 displacement, rotation : dict of csdl.Variable 78 Per-beam nodal displacements/rotations, shape (num_nodes, 3), set by ``solve()``. 79 num, dim : int 80 Number of frame nodes and of global dofs. 81 bc_indices : np.ndarray 82 Global indices of the constrained dofs. 83 """ 84 85 def __init__(self, 86 beams:List[Beam], 87 joints:Optional[List[Joint]] = None, 88 acc:Optional[csdl.Variable] = None): 89 90 if acc is not None and acc.shape != (6,): 91 raise ValueError("acc must have shape (6,)") 92 93 if not beams: 94 raise ValueError("beams must be added to the frame") 95 96 if len(set(map(id, beams))) != len(beams): 97 raise ValueError("each beam can be added to a frame only once") 98 99 if len({beam.name for beam in beams}) != len(beams): 100 raise ValueError("beam names must be unique within a frame") 101 102 self.beams = beams 103 self.joints = [] if joints is None else joints 104 self.acc = acc 105 self.displacement: Dict[str, csdl.Variable] = {} 106 self.rotation: Dict[str, csdl.Variable] = {} 107 self.U = None 108 109 # global node number of every node of every beam 110 self._nodes: Dict[Beam, np.ndarray] = self._number_nodes() 111 self.num = len(np.unique(np.concatenate(list(self._nodes.values())))) 112 self.dim = self.num * 6 113 self.bc_indices = self._bc_indices() 114 115 # global matrices and loads, with the boundary conditions applied 116 self.K, self.M = self._global_matrices() 117 self.F = self._global_loads() 118 119 120 def solve(self)->None: 121 """solve K U = F and extract the per-beam displacements and rotations""" 122 self.U = U = csdl.solve_linear(self.K, self.F) 123 124 for beam in self.beams: 125 u = self._nodal_displacements(beam) 126 self.displacement[beam.name] = u[:, :3] 127 self.rotation[beam.name] = u[:, 3:] 128 129 130 def compute_stress(self, beams:Optional[List[Beam]] = None)->Dict[str, csdl.Variable]: 131 """ 132 per-beam stresses from the cross-section stress models (call solve() first) 133 134 beams : the beams to evaluate (default: all); every one needs a 135 cross-section with a stress model 136 """ 137 if self.U is None: 138 raise RuntimeError("call solve() before compute_stress()") 139 140 stress = {} 141 for beam in self.beams if beams is None else beams: 142 # element displacements [u_a, u_b], shape (num_elements, 12) 143 u = self._nodal_displacements(beam) 144 element_displacements = csdl.concatenate([u[:-1], u[1:]], axis=1) 145 146 stress[beam.name] = beam.cs.stress(beam.element_loads(element_displacements)) 147 148 return stress 149 150 151 def compute_mass_properties(self)->Tuple[csdl.Variable, csdl.Variable]: 152 """center of gravity (3,) and total mass of the beams (extra inertial masses excluded)""" 153 mass, rmvec = 0, 0 154 for beam in self.beams: 155 mass += beam.mass 156 rmvec += beam.rmvec 157 158 cg = rmvec / mass 159 return cg, mass 160 161 162 def node_dofs(self, beam:Beam)->np.ndarray: 163 """global dof indices of every node of a beam, shape (num_nodes, 6)""" 164 return self._nodes[beam][:, None] * 6 + np.arange(6) 165 166 167 def _element_dofs(self, beam:Beam)->np.ndarray: 168 """global dof indices of every element of a beam, shape (num_elements, 12)""" 169 dofs = self.node_dofs(beam) 170 return np.hstack([dofs[:-1], dofs[1:]]) 171 172 173 def _nodal_displacements(self, beam:Beam)->csdl.Variable: 174 """rows of U for the nodes of a beam, shape (num_nodes, 6)""" 175 return self.U[self.node_dofs(beam).ravel().tolist()].reshape((beam.num_nodes, 6)) 176 177 178 def _number_nodes(self)->Dict[Beam, np.ndarray]: 179 """ 180 global node numbers: every beam node gets a provisional id, joined 181 nodes are merged with union-find (independent of the joint order) and 182 the merged nodes are numbered contiguously 183 """ 184 start, total = {}, 0 185 for beam in self.beams: 186 start[beam] = total 187 total += beam.num_nodes 188 189 parent = list(range(total)) 190 191 def find(i): 192 while parent[i] != i: 193 parent[i] = parent[parent[i]] 194 i = parent[i] 195 return i 196 197 for joint in self.joints: 198 for member in joint.members: 199 if member not in start: 200 raise ValueError(f"joint member {member.name!r} is not a beam of this frame") 201 202 # the joint's first member provides the merged node's id 203 root = find(start[joint.members[0]] + joint.nodes[0]) 204 for member, node in zip(joint.members[1:], joint.nodes[1:]): 205 parent[find(start[member] + node)] = root 206 207 # number the merged nodes in set order, which matches the numbering 208 # of earlier versions (the layout of K and U is unchanged) 209 ids = [find(i) for i in range(total)] 210 number = {node: n for n, node in enumerate(set(ids))} 211 212 return {beam: np.array([number[i] for i in ids[start[beam]:start[beam] + beam.num_nodes]]) 213 for beam in self.beams} 214 215 216 def _bc_indices(self)->np.ndarray: 217 """global indices of the constrained dofs (fixed: all 6, pinned: 3 translations)""" 218 indices = [np.empty(0, dtype=int)] 219 for beam in self.beams: 220 dofs = self.node_dofs(beam) 221 indices += [dofs[node] for node in beam.fixed_boundary_conditions] 222 indices += [dofs[node, :3] for node in beam.pinned_boundary_conditions] 223 224 return np.unique(np.concatenate(indices)) 225 226 227 def _global_matrices(self)->Tuple[csdl.Variable, csdl.Variable]: 228 """ 229 assemble the global stiffness/mass matrices from the element matrices; 230 the rows/columns of constrained dofs are never written and their 231 diagonal is set to 1 232 """ 233 dim = self.dim 234 is_bc = np.isin(np.arange(dim), self.bc_indices) 235 236 # global (row, col) of every entry of every element matrix 237 dofs = np.vstack([self._element_dofs(beam) for beam in self.beams]) 238 rows = np.repeat(dofs, 12, axis=1).ravel() 239 cols = np.tile(dofs, (1, 12)).ravel() 240 positions = np.where(is_bc[rows] | is_bc[cols], -1, rows * dim + cols) 241 242 init = np.diag(is_bc.astype(float)) 243 K = _scatter_add(init, positions, [beam.transformed_stiffness for beam in self.beams]) 244 M = _scatter_add(init, positions, [beam.transformed_mass for beam in self.beams]) 245 246 return K, M 247 248 249 def _global_loads(self)->csdl.Variable: 250 """ 251 assemble the global load vector from the nodal loads and any inertial 252 loads; loads on constrained dofs are dropped 253 """ 254 zeros = np.zeros(self.dim) 255 256 def positions(dofs): 257 dofs = np.concatenate([d.ravel() for d in dofs]) 258 return np.where(np.isin(dofs, self.bc_indices), -1, dofs) 259 260 # nodal loads 261 loaded = [beam for beam in self.beams if beam.loads is not None] 262 if loaded: 263 F = _scatter_add(zeros, positions([self.node_dofs(b) for b in loaded]), [b.loads for b in loaded]) 264 else: 265 F = csdl.Variable(value=zeros) 266 267 acc = self.acc 268 if acc is not None: 269 # structural inertial loads M @ acc, element by element: m_e @ [acc, acc] 270 acc12 = csdl.concatenate([acc, acc]) 271 inertial = [csdl.einsum(b.transformed_mass, acc12, action='ijk,k->ij') for b in self.beams] 272 F = F + _scatter_add(zeros, positions([self._element_dofs(b) for b in self.beams]), inertial) 273 274 # extra point masses are resolved as loads 275 massed = [beam for beam in self.beams if beam.extra_inertial_mass is not None] 276 if massed: 277 extra = [csdl.outer(b.extra_inertial_mass, acc) for b in massed] 278 F = F + _scatter_add(zeros, positions([self.node_dofs(b) for b in massed]), extra) 279 280 return F
A linear static frame model built from beams and joints.
Building the Frame numbers the frame nodes and assembles the global
stiffness and mass matrices and the load vector (with boundary conditions
applied); call solve() for the displacements and compute_stress()
for the stresses. The dof numbering is owned by the Frame, so a Beam can
be used in several Frames.
Parameters
- beams (list of Beam): The beams (unique objects with unique names), with loads and boundary conditions already added.
- joints (list of Joint, optional): Rigid connections between beams of this frame.
- acc (csdl.Variable, optional): Frame acceleration [ax, ay, az, alpha_x, alpha_y, alpha_z], shape (6,), applied as inertial loads M @ acc on the structure and extra masses.
Attributes
- K, M (csdl.Variable): Global stiffness and mass matrices, shape (dim, dim), dim = 6 * num.
- F (csdl.Variable): Global load vector, shape (dim,).
- U (csdl.Variable):
Global displacement vector, set by
solve(). - displacement, rotation (dict of csdl.Variable):
Per-beam nodal displacements/rotations, shape (num_nodes, 3), set by
solve(). - num, dim (int): Number of frame nodes and of global dofs.
- bc_indices (np.ndarray): Global indices of the constrained dofs.
Frame( beams: List[Beam], joints: Optional[List[Joint]] = None, acc: Optional[csdl_alpha.src.graph.variable.Variable] = None)
85 def __init__(self, 86 beams:List[Beam], 87 joints:Optional[List[Joint]] = None, 88 acc:Optional[csdl.Variable] = None): 89 90 if acc is not None and acc.shape != (6,): 91 raise ValueError("acc must have shape (6,)") 92 93 if not beams: 94 raise ValueError("beams must be added to the frame") 95 96 if len(set(map(id, beams))) != len(beams): 97 raise ValueError("each beam can be added to a frame only once") 98 99 if len({beam.name for beam in beams}) != len(beams): 100 raise ValueError("beam names must be unique within a frame") 101 102 self.beams = beams 103 self.joints = [] if joints is None else joints 104 self.acc = acc 105 self.displacement: Dict[str, csdl.Variable] = {} 106 self.rotation: Dict[str, csdl.Variable] = {} 107 self.U = None 108 109 # global node number of every node of every beam 110 self._nodes: Dict[Beam, np.ndarray] = self._number_nodes() 111 self.num = len(np.unique(np.concatenate(list(self._nodes.values())))) 112 self.dim = self.num * 6 113 self.bc_indices = self._bc_indices() 114 115 # global matrices and loads, with the boundary conditions applied 116 self.K, self.M = self._global_matrices() 117 self.F = self._global_loads()
def
solve(self) -> None:
120 def solve(self)->None: 121 """solve K U = F and extract the per-beam displacements and rotations""" 122 self.U = U = csdl.solve_linear(self.K, self.F) 123 124 for beam in self.beams: 125 u = self._nodal_displacements(beam) 126 self.displacement[beam.name] = u[:, :3] 127 self.rotation[beam.name] = u[:, 3:]
solve K U = F and extract the per-beam displacements and rotations
def
compute_stress( self, beams: Optional[List[Beam]] = None) -> Dict[str, csdl_alpha.src.graph.variable.Variable]:
130 def compute_stress(self, beams:Optional[List[Beam]] = None)->Dict[str, csdl.Variable]: 131 """ 132 per-beam stresses from the cross-section stress models (call solve() first) 133 134 beams : the beams to evaluate (default: all); every one needs a 135 cross-section with a stress model 136 """ 137 if self.U is None: 138 raise RuntimeError("call solve() before compute_stress()") 139 140 stress = {} 141 for beam in self.beams if beams is None else beams: 142 # element displacements [u_a, u_b], shape (num_elements, 12) 143 u = self._nodal_displacements(beam) 144 element_displacements = csdl.concatenate([u[:-1], u[1:]], axis=1) 145 146 stress[beam.name] = beam.cs.stress(beam.element_loads(element_displacements)) 147 148 return stress
per-beam stresses from the cross-section stress models (call solve() first)
beams : the beams to evaluate (default: all); every one needs a cross-section with a stress model
def
compute_mass_properties( self) -> Tuple[csdl_alpha.src.graph.variable.Variable, csdl_alpha.src.graph.variable.Variable]:
151 def compute_mass_properties(self)->Tuple[csdl.Variable, csdl.Variable]: 152 """center of gravity (3,) and total mass of the beams (extra inertial masses excluded)""" 153 mass, rmvec = 0, 0 154 for beam in self.beams: 155 mass += beam.mass 156 rmvec += beam.rmvec 157 158 cg = rmvec / mass 159 return cg, mass
center of gravity (3,) and total mass of the beams (extra inertial masses excluded)
def
mesh_from_points_and_edges( points: numpy.ndarray, edges: numpy.ndarray, num_nodes: int) -> numpy.ndarray:
8def mesh_from_points_and_edges(points:np.ndarray, edges:np.ndarray, num_nodes:int)->np.ndarray: 9 """ 10 Straight beam meshes along the edges of a point cloud. 11 12 Parameters 13 ---------- 14 points : np.ndarray 15 Point coordinates, shape (num_points, 3). 16 edges : np.ndarray 17 Point index pairs (start, end), shape (num_edges, 2). 18 num_nodes : int 19 Number of evenly spaced nodes per edge (including both end points). 20 21 Returns 22 ------- 23 np.ndarray 24 Node coordinates, shape (num_edges, num_nodes, 3). 25 """ 26 points = np.asarray(points) 27 edges = np.asarray(edges) 28 return np.linspace(points[edges[:, 0]], points[edges[:, 1]], num=num_nodes, axis=1)
Straight beam meshes along the edges of a point cloud.
Parameters
- points (np.ndarray): Point coordinates, shape (num_points, 3).
- edges (np.ndarray): Point index pairs (start, end), shape (num_edges, 2).
- num_nodes (int): Number of evenly spaced nodes per edge (including both end points).
Returns
- np.ndarray: Node coordinates, shape (num_edges, num_nodes, 3).
def
plot_box( plotter, mesh, height, width, cell_data=None, cmap='viridis', color='lightblue'):
32def plot_box(plotter, mesh, height, width, cell_data=None, cmap='viridis', color='lightblue'): 33 """ 34 draw a rectangular box per element 35 36 height, width : per-element box dimensions 37 cell_data : optional per-element values used for coloring 38 """ 39 n = mesh.shape[0] 40 colors = _colorize(cell_data, cmap, color, n) 41 42 up = np.array([0, 0, 1]) 43 for i in range(n - 1): 44 start = mesh[i, :] 45 end = mesh[i + 1, :] 46 47 direction = end - start 48 length = np.linalg.norm(direction) 49 direction /= length 50 51 # rotate the z-aligned cube onto the element axis 52 if np.allclose(direction, up): 53 axis = up 54 angle = 0 55 else: 56 axis = np.cross(up, direction) 57 angle = np.rad2deg(np.arccos(np.dot(up, direction))) 58 59 midpoint = (start + end) / 2 60 box = pv.Cube(center=midpoint, x_length=width[i], y_length=height[i], z_length=length) 61 box.rotate_vector(vector=axis, angle=angle, point=midpoint, inplace=True) 62 63 plotter.add_mesh(box, color=colors[i], show_edges=True)
draw a rectangular box per element
height, width : per-element box dimensions cell_data : optional per-element values used for coloring
def
plot_cyl( plotter, mesh, radius, cell_data=None, cmap='viridis', color='lightblue'):
66def plot_cyl(plotter, mesh, radius, cell_data=None, cmap='viridis', color='lightblue'): 67 """ 68 draw a cylinder per element 69 70 radius : per-element cylinder radius 71 cell_data : optional per-element values used for coloring 72 """ 73 n = mesh.shape[0] 74 colors = _colorize(cell_data, cmap, color, n) 75 76 for i in range(n - 1): 77 start = mesh[i, :] 78 end = mesh[i + 1, :] 79 80 cyl = pv.Cylinder(center=(start + end) / 2, 81 direction=end - start, 82 radius=radius[i], 83 height=np.linalg.norm(end - start)) 84 85 plotter.add_mesh(cyl, color=colors[i])
draw a cylinder per element
radius : per-element cylinder radius cell_data : optional per-element values used for coloring
def
plot_mesh( plotter, mesh, cell_data=None, plot_mesh=True, cmap='viridis', line_width=20, render_lines_as_tubes=True, color='red'):
88def plot_mesh(plotter, 89 mesh, 90 cell_data=None, 91 plot_mesh=True, 92 cmap='viridis', 93 line_width=20, 94 render_lines_as_tubes=True, 95 color='red'): 96 """ 97 draw the beam as a polyline, optionally colored by per-element cell_data 98 (plot_mesh is unused and kept for backward compatibility) 99 """ 100 polyline = _polyline(mesh) 101 102 if cell_data is not None: 103 # normalize the cell data 104 polyline.cell_data['cell_data'] = (cell_data - cell_data.min()) / (cell_data.max() - cell_data.min()) 105 plotter.add_mesh(polyline, 106 scalars='cell_data', 107 render_lines_as_tubes=render_lines_as_tubes, 108 style='wireframe', 109 line_width=line_width, 110 cmap=cmap, 111 show_scalar_bar=False) 112 else: 113 plotter.add_mesh(polyline, 114 render_lines_as_tubes=render_lines_as_tubes, 115 style='wireframe', 116 line_width=line_width, 117 color=color, 118 show_scalar_bar=False)
draw the beam as a polyline, optionally colored by per-element cell_data (plot_mesh is unused and kept for backward compatibility)
def
plot_points( plotter, mesh, color='red', point_size=50, render_points_as_spheres=True):
121def plot_points(plotter, mesh, color='red', point_size=50, render_points_as_spheres=True): 122 """draw the mesh nodes as points""" 123 plotter.add_points(np.asarray(mesh, dtype=float), 124 color=color, 125 point_size=point_size, 126 render_points_as_spheres=render_points_as_spheres)
draw the mesh nodes as points