| 39 | |
| 40 | @dataclass |
| 41 | class Rod_Elem: |
| 42 | id : str |
| 43 | node_a: Node |
| 44 | node_b: Node |
| 45 | radius: float |
| 46 | E : float # Young's modulus of elasticity |
| 47 | rho : float # density |
| 48 | # Internal values computed in __post_init__() |
| 49 | # These are not supplied by the caller when |
| 50 | # creating a new Rod_Elem; instead they are |
| 51 | # computed from the above inputs. |
| 52 | dX : float = field(init=False) |
| 53 | dY : float = field(init=False) |
| 54 | length : float = field(init=False) |
| 55 | cross_sect_area: float = field(init=False) |
| 56 | mass : float = field(init=False) |
| 57 | |
| 58 | def __post_init__(self): |
| 59 | """ |
| 60 | Called when a new Rod_Elem is defined via |
| 61 | element = Rod_Elem(id, node_a, node_b, E, rho) |
| 62 | """ |
| 63 | self.dX = self.node_b.x - self.node_a.x |
| 64 | self.dY = self.node_b.y - self.node_a.y |
| 65 | self.length = np.sqrt( self.dX**2 + self.dY**2 ) |
| 66 | self.cross_sect_area = np.pi * self.radius**2 |
| 67 | self.mass = self.rho * self.length * self.cross_sect_area |
| 68 | |
| 69 | def stiffness_matrix(self): |
| 70 | """ |
| 71 | http://www.ita.uni-heidelberg.de/~dullemond/lectures/num_phys_2010/Chapter_FiniteElements.pdf |
| 72 | """ |
| 73 | K = self.E * self.cross_sect_area / self.length |
| 74 | cos = self.dX/self.length |
| 75 | sin = self.dY/self.length |
| 76 | cos2 = cos**2 |
| 77 | sin2 = sin**2 |
| 78 | sincos = sin*cos |
| 79 | return K*np.array([[ cos2 , sincos, -cos2 , -sincos ], |
| 80 | [ sincos, sin2 , -sincos, -sin2 ], |
| 81 | [-cos2 , -sincos, cos2 , sincos, ], |
| 82 | [-sincos, -sin2 , sincos, sin2 , ],]) |
| 83 | |
| 84 | def mass_matrix(self): |
| 85 | """ |
| 86 | Lumped mass, returned as the diagonal of the matrix. |
| 87 | http://kis.tu.kielce.pl/mo/COLORADO_FEM/colorado/IFEM.Ch31.pdf, p 31-7 |
| 88 | """ |
| 89 | return 0.5 * self.mass * np.array([1.0, 1.0, 1.0, 1.0]) |
| 90 | |
| 91 | def main(): |
| 92 | p_a = Node('A', 1.0, 1.0) |