1
+ {
2
+ "nbformat" : 4 ,
3
+ "nbformat_minor" : 0 ,
4
+ "metadata" : {
5
+ "colab" : {
6
+ "private_outputs" : true ,
7
+ "provenance" : [],
8
+ "authorship_tag" : " ABX9TyNYWx+I2cFrt4q7k9DFe8h8" ,
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/1DElement.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" : " EVc6sM8n_6WC"
35
+ },
36
+ "outputs" : [],
37
+ "source" : [
38
+ " import numpy as np\n " ,
39
+ " import matplotlib.pyplot as plt\n " ,
40
+ " \n " ,
41
+ " # Define the function to assemble the global stiffness matrix\n " ,
42
+ " def assemble_stiffness_matrix(nodal_coords, nodal_connectivity, material_props):\n " ,
43
+ " num_nodes = len(nodal_coords)\n " ,
44
+ " K_global = np.zeros((num_nodes, num_nodes))\n " ,
45
+ " \n " ,
46
+ " for element in nodal_connectivity:\n " ,
47
+ " node1, node2 = element\n " ,
48
+ " x1, x2 = nodal_coords[node1], nodal_coords[node2]\n " ,
49
+ " length = x2 - x1\n " ,
50
+ " k = material_props['E'] * material_props['A'] / length\n " ,
51
+ " \n " ,
52
+ " # Element stiffness matrix for 1D linear element\n " ,
53
+ " K_local = (k / length) * np.array([[1, -1], [-1, 1]])\n " ,
54
+ " \n " ,
55
+ " # Assembly into the global stiffness matrix\n " ,
56
+ " K_global[node1:node2+1, node1:node2+1] += K_local\n " ,
57
+ " \n " ,
58
+ " return K_global\n " ,
59
+ " \n " ,
60
+ " # Apply boundary conditions (Dirichlet)\n " ,
61
+ " def apply_boundary_conditions(K, F, bc):\n " ,
62
+ " for node, value in bc.items():\n " ,
63
+ " # Modify stiffness matrix and force vector to enforce boundary condition\n " ,
64
+ " K[node, :] = 0\n " ,
65
+ " K[:, node] = 0\n " ,
66
+ " K[node, node] = 1\n " ,
67
+ " F[node] = value\n " ,
68
+ " return K, F\n " ,
69
+ " \n " ,
70
+ " # Define function to solve FEM problem\n " ,
71
+ " def fem_1d(nodal_coords, nodal_connectivity, material_props, boundary_conditions, load):\n " ,
72
+ " num_nodes = len(nodal_coords)\n " ,
73
+ " \n " ,
74
+ " # Step 1: Assemble global stiffness matrix\n " ,
75
+ " K_global = assemble_stiffness_matrix(nodal_coords, nodal_connectivity, material_props)\n " ,
76
+ " \n " ,
77
+ " # Step 2: Assemble global force vector\n " ,
78
+ " F_global = np.zeros(num_nodes)\n " ,
79
+ " for node, value in load.items():\n " ,
80
+ " F_global[node] = value\n " ,
81
+ " \n " ,
82
+ " # Step 3: Apply boundary conditions\n " ,
83
+ " K_global, F_global = apply_boundary_conditions(K_global, F_global, boundary_conditions)\n " ,
84
+ " \n " ,
85
+ " # Step 4: Solve the system of equations\n " ,
86
+ " displacements = np.linalg.solve(K_global, F_global)\n " ,
87
+ " \n " ,
88
+ " return displacements\n " ,
89
+ " \n " ,
90
+ " # Example Input Data\n " ,
91
+ " nodal_coords = np.array([0.0, 0.5, 1.0]) # Coordinates of nodes\n " ,
92
+ " nodal_connectivity = np.array([[0, 1], [1, 2]]) # Elements defined by node numbers\n " ,
93
+ " material_props = {'E': 210e9, 'A': 1e-4} # Material properties: E (Young's modulus) and A (Cross-sectional area)\n " ,
94
+ " boundary_conditions = {0: 0.0} # Boundary condition: node 0 is fixed (Dirichlet condition)\n " ,
95
+ " load = {2: 1000.0} # Load applied at node 2\n " ,
96
+ " \n " ,
97
+ " # Solve the FEM problem\n " ,
98
+ " displacements = fem_1d(nodal_coords, nodal_connectivity, material_props, boundary_conditions, load)\n " ,
99
+ " \n " ,
100
+ " # Output displacements\n " ,
101
+ " print(\" Nodal Displacements:\" , displacements)\n " ,
102
+ " \n " ,
103
+ " # Plot the displacements\n " ,
104
+ " plt.plot(nodal_coords, displacements, '-o')\n " ,
105
+ " plt.xlabel('Position (m)')\n " ,
106
+ " plt.ylabel('Displacement (m)')\n " ,
107
+ " plt.title('1D FEM Nodal Displacements')\n " ,
108
+ " plt.grid(True)\n " ,
109
+ " plt.show()\n "
110
+ ]
111
+ }
112
+ ]
113
+ }
0 commit comments