1
+ {
2
+ "nbformat" : 4 ,
3
+ "nbformat_minor" : 0 ,
4
+ "metadata" : {
5
+ "colab" : {
6
+ "private_outputs" : true ,
7
+ "provenance" : [],
8
+ "authorship_tag" : " ABX9TyMxahuOixvTTe6EdHJBA9EE" ,
9
+ "include_colab_link" : true
10
+ },
11
+ "kernelspec" : {
12
+ "name" : " python3" ,
13
+ "display_name" : " Python 3"
14
+ },
15
+ "language_info" : {
16
+ "name" : " python"
17
+ }
18
+ },
19
+ "cells" : [
20
+ {
21
+ "cell_type" : " markdown" ,
22
+ "metadata" : {
23
+ "id" : " view-in-github" ,
24
+ "colab_type" : " text"
25
+ },
26
+ "source" : [
27
+ " <a href=\" https://colab.research.google.com/github/GEORMC/Nnumerical_Methods_Course/blob/main/Truss_Element.ipynb\" target=\" _parent\" ><img src=\" https://colab.research.google.com/assets/colab-badge.svg\" alt=\" Open In Colab\" /></a>"
28
+ ]
29
+ },
30
+ {
31
+ "cell_type" : " code" ,
32
+ "execution_count" : null ,
33
+ "metadata" : {
34
+ "id" : " 7_YgV2riW2t1"
35
+ },
36
+ "outputs" : [],
37
+ "source" : [
38
+ " import numpy as np\n " ,
39
+ " # Input data\n " ,
40
+ " coordinates = np.array([\n " ,
41
+ " [0, 0],\n " ,
42
+ " [4, 0],\n " ,
43
+ " [8, 0],\n " ,
44
+ " [4, -6]\n " ,
45
+ " ])\n " ,
46
+ " connectivity = np.array([\n " ,
47
+ " [0, 3],\n " ,
48
+ " [1, 3],\n " ,
49
+ " [2, 3]\n " ,
50
+ " ])\n " ,
51
+ " E = 1 # Young's Modulus (Pa)\n " ,
52
+ " A = 1 # Cross-sectional area (m^2)\n " ,
53
+ " # Define supports (0 for not supported, 1 for supported)\n " ,
54
+ " supports = np.array([\n " ,
55
+ " [1, 1],\n " ,
56
+ " [1, 1],\n " ,
57
+ " [1, 1],\n " ,
58
+ " [0, 0]\n " ,
59
+ " ])\n " ,
60
+ " # Define applied loads (0 for no load, specify direction and value)\n " ,
61
+ " applied_loads = np.array([\n " ,
62
+ " [0, 0],\n " ,
63
+ " [0, 0], # Node 1, applied horizontal load\n " ,
64
+ " [0, 0], # Node 2, applied vertical load\n " ,
65
+ " [100, -100]\n " ,
66
+ " ])\n " ,
67
+ " # Create freedom matrix\n " ,
68
+ " num_nodes = len(coordinates)\n " ,
69
+ " num_dofs = 2 * num_nodes\n " ,
70
+ " num_elements=len(connectivity)\n " ,
71
+ " dofs = np.zeros((num_elements, 4), dtype=int)\n " ,
72
+ " NodeDof=np.zeros((num_nodes, 2) , dtype=int)\n " ,
73
+ " # Initialize global stiffness matrix and force vector\n " ,
74
+ " KG = np.zeros((num_dofs, num_dofs))\n " ,
75
+ " kt_global = np.zeros((num_dofs, num_dofs))\n " ,
76
+ " F_global = np.zeros(num_dofs)\n " ,
77
+ " # Calculate Dof Matrix\n " ,
78
+ " for i, (node1, node2) in enumerate(connectivity):\n " ,
79
+ " dofs[i,:] = np.array([2 * node1, 2 * node1 + 1, 2 * node2, 2 * node2 + 1])\n " ,
80
+ " NodeDof[node1,0]=np.array([2 * node1])\n " ,
81
+ " NodeDof[node1,1]=np.array([2 * node1 + 1])\n " ,
82
+ " NodeDof[node2,0]=np.array([ 2 * node2])\n " ,
83
+ " NodeDof[node2,1]=np.array([2* node2+1])\n " ,
84
+ " print(NodeDof)\n " ,
85
+ " # Calculate element lengths and stiffness matrices\n " ,
86
+ " for i, (node1, node2) in enumerate(connectivity):\n " ,
87
+ " x1, y1 = coordinates[node1]\n " ,
88
+ " x2, y2 = coordinates[node2]\n " ,
89
+ " L = np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)\n " ,
90
+ " c = (x2 - x1) / L\n " ,
91
+ " s = (y2 - y1) / L\n " ,
92
+ " k_local = (E * A / L) * np.array([\n " ,
93
+ " [c ** 2, c * s, -c ** 2, -c * s],\n " ,
94
+ " [c * s, s ** 2, -c * s, -s ** 2],\n " ,
95
+ " [-c ** 2, -c * s, c ** 2, c * s],\n " ,
96
+ " [-c * s, -s ** 2, c * s, s ** 2]\n " ,
97
+ " ])\n " ,
98
+ " # Assemble local stiffness matrix into global stiffness matrix\n " ,
99
+ " \n " ,
100
+ " KG[np.ix_(dofs[i,:], dofs[i,:])] += k_local\n " ,
101
+ " # print(k_local)\n " ,
102
+ " # print(np.ix_(dofs[i,:], dofs[i,:]) )\n " ,
103
+ " # print(K_global)\n " ,
104
+ " print(KG)\n " ,
105
+ " kt_global=KG*1\n " ,
106
+ " # Apply boundary conditions and applied loads\n " ,
107
+ " for node in range(num_nodes):\n " ,
108
+ " for i in range(2):\n " ,
109
+ " if supports[node, i] == 1:\n " ,
110
+ " fixed_dof = NodeDof[node,i]\n " ,
111
+ " kt_global[fixed_dof, :] = 0\n " ,
112
+ " kt_global[:, fixed_dof] = 0\n " ,
113
+ " kt_global[fixed_dof, fixed_dof] = 1\n " ,
114
+ " F_global[fixed_dof] = 0\n " ,
115
+ " else:\n " ,
116
+ " applied_force = applied_loads[node, i]\n " ,
117
+ " F_global[2 * node + i] = applied_force\n " ,
118
+ " print(KG)\n " ,
119
+ " # Solve for displacements\n " ,
120
+ " displacement = np.linalg.solve(kt_global, F_global)\n " ,
121
+ " print(KG)\n " ,
122
+ " # Calculate reaction forces\n " ,
123
+ " reaction_forces = np.dot(KG, displacement)\n " ,
124
+ " print(KG)\n " ,
125
+ " print(\" Displacements (mm):\" )\n " ,
126
+ " print(displacement * 1000)\n " ,
127
+ " print(\" Reaction Forces (N):\" )\n " ,
128
+ " print(reaction_forces)\n " ,
129
+ " \n " ,
130
+ " \n "
131
+ ]
132
+ }
133
+ ]
134
+ }
0 commit comments