# ***********************************
# Author: Pedro Jorge De Los Santos
# E-mail: delossantosmfq@gmail.com
# Blog: numython.github.io
# License: MIT License
# ***********************************
import numpy as np
from .node import Node
[docs]
class Element:
"""Base class for finite elements.
Elements own formulation, connectivity, and physical properties. Solved
response belongs to analysis results, not to the element instance.
"""
def __init__(self, etype):
self.etype = etype
self.label = None
def __str__(self):
return str(self.__class__)
#~ =========================== NODE ===========================
def _validate_nodes(nodes, expected_count, element_name):
"""Return validated element connectivity as a tuple of finite 2D nodes."""
try:
nodes = tuple(nodes)
except TypeError as exc:
raise ValueError(
f"{element_name} nodes must be an iterable of Node objects"
) from exc
if len(nodes) != expected_count:
raise ValueError(
f"{element_name} requires exactly {expected_count} nodes; "
f"got {len(nodes)}"
)
if not all(isinstance(node, Node) for node in nodes):
raise ValueError(f"{element_name} connectivity must contain only Node objects")
if not all(np.isfinite([node.x, node.y]).all() for node in nodes):
raise ValueError(f"{element_name} node coordinates must be finite")
return nodes
def _positive_finite(value, name):
"""Return a finite positive scalar material or section property."""
try:
value = float(value)
except (TypeError, ValueError) as exc:
raise ValueError(f"{name} must be a finite positive scalar") from exc
if not np.isfinite(value) or value <= 0.0:
raise ValueError(f"{name} must be a finite positive scalar")
return value
def _poisson_ratio(value):
"""Return a physically admissible isotropic Poisson ratio."""
try:
value = float(value)
except (TypeError, ValueError) as exc:
raise ValueError("nu must be a finite scalar in the range -1 < nu < 0.5") from exc
if not np.isfinite(value) or not (-1.0 < value < 0.5):
raise ValueError("nu must be a finite scalar in the range -1 < nu < 0.5")
return value
def _validate_nonzero_length(nodes, element_name):
"""Reject coincident-node connectivity for length-dependent elements."""
n1, n2 = nodes
if np.hypot(n2.x - n1.x, n2.y - n1.y) == 0.0:
raise ValueError(f"{element_name} requires two distinct node coordinates")
def _validate_triangle_geometry(nodes):
"""Reject collinear connectivity for the constant-strain triangle."""
n1, n2, n3 = nodes
twice_area = (
n1.x * (n2.y - n3.y)
+ n2.x * (n3.y - n1.y)
+ n3.x * (n1.y - n2.y)
)
if twice_area == 0.0:
raise ValueError("LinearTriangle requires three non-collinear nodes")
[docs]
class Spring(Element):
"""
Spring element for finite element analysis
*nodes* : tuple
Connectivity for element given as tuple of
:class:`~nusa.node.Node` objects
*ke* : float
Spring stiffness
Example ::
n1 = Node((0,0))
n2 = Node((0,0))
e1 = Spring((n1,n2), 1000)
"""
def __init__(self,nodes,ke):
Element.__init__(self, etype="spring")
self.nodes = _validate_nodes(nodes, 2, "Spring")
self.k = _positive_finite(ke, "Spring stiffness k")
[docs]
def compute_results(self, u_e):
"""Return canonical spring results from local displacements."""
u_e = np.asarray(u_e, dtype=float).reshape(-1)
if u_e.size != 2:
raise ValueError("Spring result evaluation requires 2 displacements")
values = self.get_element_stiffness() @ u_e
return {
"force_i": float(values[0]),
"force_j": float(values[1]),
}
[docs]
def get_element_stiffness(self):
r"""
Get stiffness matrix for this element.
The stiffness matrix for a spring element is defined by:
.. math::
[k]_e = \begin{bmatrix}
k & -k \\
-k & k \\
\end{bmatrix}
where *k* is the spring stiffness.
Return a numpy array.
"""
self._KE = np.array([[self.k,-self.k],[-self.k,self.k]])
return self._KE
[docs]
class Bar(Element):
"""
Bar element for finite element analysis
*nodes* : :class:`~nusa.node.Node`
Connectivity for element
*E* : float
Young's modulus
*A* : float
Area of element
"""
def __init__(self,nodes,E,A):
Element.__init__(self,etype="bar")
self.nodes = _validate_nodes(nodes, 2, "Bar")
_validate_nonzero_length(self.nodes, "Bar")
self.E = _positive_finite(E, "Young's modulus E")
self.A = _positive_finite(A, "Cross-sectional area A")
[docs]
def compute_results(self, u_e):
"""Return canonical bar results from local displacements."""
u_e = np.asarray(u_e, dtype=float).reshape(-1)
if u_e.size != 2:
raise ValueError("Bar result evaluation requires 2 displacements")
forces = self.get_element_stiffness() @ u_e
axial_strain = (u_e[1] - u_e[0]) / self.L
axial_force = self.E * self.A * axial_strain
return {
"force_i": float(forces[0]),
"force_j": float(forces[1]),
"axial_force": float(axial_force),
"axial_stress": float(axial_force / self.A),
}
@property
def L(self):
"""
Length of element
"""
ni,nj = self.nodes
x0,x1,y0,y1 = ni.x, nj.x, ni.y, nj.y
_l = np.sqrt( (x1-x0)**2 + (y1-y0)**2 )
return _l
[docs]
def get_element_stiffness(self):
r"""
Get stiffness matrix for this element
The stiffness matrix for bar element is given by:
.. math::
[k]_e = \frac{AE}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}
where
* A - Cross-section of element
* E - Young's Modulus
* L - Length of element
"""
self._KE = (self.A*self.E/self.L)*np.array([[1,-1],[-1,1]])
return self._KE
[docs]
class Truss(Element):
"""
Truss element for finite element analysis
*nodes* : Tuple of :class:`~nusa.node.Node`
Connectivity for element
*E* : float
Young modulus
*A* : float
Area of element
"""
def __init__(self,nodes,E,A):
Element.__init__(self,etype="truss")
self.nodes = _validate_nodes(nodes, 2, "Truss")
_validate_nonzero_length(self.nodes, "Truss")
self.E = _positive_finite(E, "Young's modulus E")
self.A = _positive_finite(A, "Cross-sectional area A")
@property
def L(self):
"""
Length of element
"""
ni,nj = self.nodes
x0,x1,y0,y1 = ni.x, nj.x, ni.y, nj.y
_l = np.sqrt( (x1-x0)**2 + (y1-y0)**2 )
return _l
@property
def theta(self):
"""
Element angle, measure from X-positive axis counter-clockwise.
"""
ni,nj = self.nodes
x0,x1,y0,y1 = ni.x, nj.x, ni.y, nj.y
# ~ if x0==x1:
# ~ theta = 90*(np.pi/180)
# ~ else:
theta = np.arctan2((y1-y0),(x1-x0))
return theta
[docs]
def compute_results(self, u_e):
"""Return canonical truss results from global element displacements."""
u_e = np.asarray(u_e, dtype=float).reshape(-1)
if u_e.size != 4:
raise ValueError("Truss result evaluation requires 4 displacements")
C = np.cos(self.theta)
S = np.sin(self.theta)
axial_force = (self.E * self.A / self.L) * np.dot(
np.array([-C, -S, C, S], dtype=float),
u_e,
)
return {
"axial_force": float(axial_force),
"axial_stress": float(axial_force / self.A),
}
[docs]
def get_element_stiffness(self):
"""
Get stiffness matrix for this element
"""
multiplier = (self.A*self.E/self.L)
C = np.cos(self.theta)
S = np.sin(self.theta)
CS = C*S
self._K = multiplier*np.array([[C**2 , CS , -C**2, -CS ],
[CS , S**2 , -CS , -S**2],
[-C**2, -CS , C**2 , CS ],
[-CS , -S**2, CS , S**2 ]])
return self._K
[docs]
class Beam(Element):
"""
Beam element for finite element analysis
*nodes* : :class:`~nusa.node.Node`
Connectivity for element
*E* : float
Young's modulus
*I* : float
Moment of inertia
"""
def __init__(self,nodes,E,I):
Element.__init__(self,etype="beam")
self.nodes = _validate_nodes(nodes, 2, "Beam")
_validate_nonzero_length(self.nodes, "Beam")
self.E = _positive_finite(E, "Young's modulus E")
self.I = _positive_finite(I, "Second moment of area I")
[docs]
def get_element_stiffness(self):
"""
Get stiffness matrix for this element
"""
multiplier = (self.I*self.E/self.L**3)
a = 6*self.L
b = 4*self.L**2
c = 2*self.L**2
self._K = multiplier*np.array([[ 12, a, -12, a],
[ a, b, -a, c],
[-12,-a, 12,-a],
[ a, c, -a, b]])
return self._K
[docs]
def compute_results(self, u_e):
"""Return canonical beam end actions from local displacements."""
u_e = np.asarray(u_e, dtype=float).reshape(-1)
if u_e.size != 4:
raise ValueError("Beam result evaluation requires 4 displacements")
actions = self.get_element_stiffness() @ u_e
return {
"shear_force_i": float(actions[0]),
"shear_force_j": float(actions[2]),
"bending_moment_i": float(actions[1]),
"bending_moment_j": float(actions[3]),
}
@property
def L(self):
"""
Length of element
"""
ni,nj = self.nodes
x0,x1,y0,y1 = ni.x, nj.x, ni.y, nj.y
_l = np.sqrt( (x1-x0)**2 + (y1-y0)**2 )
return _l
[docs]
class LinearTriangle(Element):
"""
Linear triangle element for finite element analysis
*nodes* : :class:`~nusa.node.Node`
Connectivity for element
*E* : float
Young's modulus
*nu* : float
Poisson ratio
*t* : float
Thickness
Example::
n1 = Node((0,0))
n2 = Node((0.5,0))
n3 = Node((0.5,0.25))
e1 = LinearTriangle((n1,n2,n3),210e9, 0.3, 0.025)
"""
def __init__(self,nodes,E,nu,t):
Element.__init__(self,etype="triangle")
self.nodes = _validate_nodes(nodes, 3, "LinearTriangle")
_validate_triangle_geometry(self.nodes)
self.E = _positive_finite(E, "Young's modulus E")
self.nu = _poisson_ratio(nu)
self.t = _positive_finite(t, "Thickness t")
@property
def D(self):
"""
Constitutive matrix
Currently only plane stress supported
"""
nu, E = self.nu, self.E
D = (E/(1-nu**2))*np.array([[1, nu, 0],
[nu, 1, 0],
[0, 0, (1-nu)/2]
])
return D
def _signed_area(self):
n1, n2, n3 = self.nodes
xi, yi = n1.x, n1.y
xj, yj = n2.x, n2.y
xm, ym = n3.x, n3.y
return (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj))/2
@property
def B(self):
ni, nj, nm = self.nodes
signed_area = self._signed_area()
if signed_area == 0.0:
raise ValueError("Degenerate triangle: area must be nonzero")
betai = nj.y - nm.y
betaj = nm.y - ni.y
betam = ni.y - nj.y
gammai = nm.x - nj.x
gammaj = ni.x - nm.x
gammam = nj.x - ni.x
B = (1/(2*signed_area))*np.array([[betai, 0, betaj, 0, betam, 0],
[0, gammai, 0, gammaj, 0, gammam],
[gammai, betai, gammaj, betaj, gammam, betam]
])
return B
@property
def A(self):
return abs(self._signed_area())
[docs]
def get_element_stiffness(self):
"""
Get stiffness matrix for this element
"""
ni, nj, nm = self.nodes
A, nu, t, E = self.A, self.nu, self.t, self.E
B, D = self.B, self.D
return t*A*np.dot(np.dot(B.T,D),B)
[docs]
def compute_strain(self, u_e):
"""Return the constant engineering strain vector for local displacements."""
u_e = np.asarray(u_e, dtype=float).reshape(-1)
if u_e.size != 6:
raise ValueError(
"LinearTriangle strain evaluation requires 6 displacements"
)
return self.B @ u_e
[docs]
def compute_stress(self, u_e):
"""Return the constant stress vector for local displacements."""
return self.D @ self.compute_strain(u_e)
[docs]
def compute_results(self, u_e):
"""Return canonical CST stress/strain results from local displacements."""
strain = np.asarray(self.compute_strain(u_e), dtype=float).reshape(-1)
stress = np.asarray(self.D @ strain, dtype=float).reshape(-1)
return {
"stress_xx": float(stress[0]),
"stress_yy": float(stress[1]),
"stress_xy": float(stress[2]),
"strain_xx": float(strain[0]),
"strain_yy": float(strain[1]),
"strain_xy": float(strain[2]),
}
if __name__=='__main__':
pass