# ***********************************
# Author: Pedro Jorge De Los Santos
# E-mail: delossantosmfq@gmail.com
# Blog: numython.github.io
# License: MIT License
# ***********************************
import numpy as np
from .material import Material
from .node import Node
from .section import Section
[docs]
class Element:
"""Base class for finite elements.
Elements own formulation and connectivity and reference the physical
properties required by that formulation. 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."""
if isinstance(value, (bool, np.bool_)):
raise ValueError(f"{name} must be a finite positive scalar")
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 _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.
Parameters
----------
nodes : tuple of Node
Element connectivity.
material : Material
Material providing Young's modulus E.
section : Section
Section providing cross-sectional area A.
'''
def __init__(self, nodes, *, material, section):
Element.__init__(self, etype='bar')
self.nodes = _validate_nodes(nodes, 2, 'Bar')
_validate_nonzero_length(self.nodes, 'Bar')
if not isinstance(material, Material):
raise TypeError('Bar material must be a Material instance')
if not isinstance(section, Section):
raise TypeError('Bar section must be a Section instance')
if section.A is None:
raise ValueError("Bar requires section property 'A'")
self._material = material
self._section = section
@property
def material(self):
"""Material associated with the element (read-only)."""
return self._material
@property
def section(self):
"""Section associated with the element (read-only)."""
return self._section
@property
def E(self):
'''Young's modulus provided by the material.'''
return self.material.E
@property
def A(self):
'''Cross-sectional area provided by the section.'''
return self.section.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
return np.sqrt((x1-x0)**2 + (y1-y0)**2)
[docs]
def get_element_stiffness(self):
'''Return the bar element stiffness matrix.'''
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.
Parameters
----------
nodes : tuple of Node
Element connectivity.
material : Material
Material providing Young's modulus E.
section : Section
Section providing cross-sectional area A.
'''
def __init__(self, nodes, *, material, section):
Element.__init__(self, etype='truss')
self.nodes = _validate_nodes(nodes, 2, 'Truss')
_validate_nonzero_length(self.nodes, 'Truss')
if not isinstance(material, Material):
raise TypeError('Truss material must be a Material instance')
if not isinstance(section, Section):
raise TypeError('Truss section must be a Section instance')
if section.A is None:
raise ValueError("Truss requires section property 'A'")
self._material = material
self._section = section
@property
def material(self):
"""Material associated with the element (read-only)."""
return self._material
@property
def section(self):
"""Section associated with the element (read-only)."""
return self._section
@property
def E(self):
'''Young's modulus provided by the material.'''
return self.material.E
@property
def A(self):
'''Cross-sectional area provided by the section.'''
return self.section.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):
'''Euler-Bernoulli beam element for finite element analysis.
Parameters
----------
nodes : tuple of Node
Element connectivity.
material : Material
Material providing Young's modulus E.
section : Section
Section providing second moment of area I.
'''
def __init__(self, nodes, *, material, section):
Element.__init__(self, etype='beam')
self.nodes = _validate_nodes(nodes, 2, 'Beam')
_validate_nonzero_length(self.nodes, 'Beam')
if not isinstance(material, Material):
raise TypeError('Beam material must be a Material instance')
if not isinstance(section, Section):
raise TypeError('Beam section must be a Section instance')
if section.I is None:
raise ValueError("Beam requires section property 'I'")
self._material = material
self._section = section
@property
def material(self):
"""Material associated with the element (read-only)."""
return self._material
@property
def section(self):
"""Section associated with the element (read-only)."""
return self._section
@property
def E(self):
'''Young's modulus provided by the material.'''
return self.material.E
@property
def I(self):
'''Second moment of area provided by the section.'''
return self.section.I
[docs]
def get_element_stiffness(self):
'''Return the beam element stiffness matrix.'''
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
return np.sqrt((x1-x0)**2 + (y1-y0)**2)
[docs]
class LinearTriangle(Element):
'''Three-node constant-strain triangle (plane stress).
Parameters
----------
nodes : tuple of Node
Triangle connectivity.
material : Material
Material defining E and nu (both required).
thickness : float
Finite, strictly positive thickness.
'''
def __init__(self, nodes, *, material, thickness):
Element.__init__(self, etype='triangle')
self.nodes = _validate_nodes(nodes, 3, 'LinearTriangle')
_validate_triangle_geometry(self.nodes)
if not isinstance(material, Material):
raise TypeError('LinearTriangle material must be a Material instance')
if material.nu is None:
raise ValueError("LinearTriangle requires material property 'nu'")
self._material = material
self._thickness = _positive_finite(thickness, 'Thickness')
@property
def material(self):
return self._material
@property
def E(self):
return self.material.E
@property
def nu(self):
return self.material.nu
@property
def thickness(self):
return self._thickness
@property
def t(self):
return self.thickness
@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