'Main'-funksjon
I main har vi laget en hovedramme for hvordan programmet kan bygges opp. Her er også forslag til flere av funksjonskallene inkludert slik at det blir lettere å komme igang med programmeringen. Legg merke til at funksjonene lesinput, lengder og all visualisering allerede er implementert. Flere detaljer om disse funksjonene kommer i neste avsnitt.
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
from structure_visualization import * # -----Rammeanalyse----- def main(): # -----Initialiserer figurer----- fig_init, ax_init, fig_def, ax_def = setup_plots() # -----Til visualiseringen, velg første indeks brukt i nummerering av noder og element----- first_index = 1 # -----Leser input-data----- npunkt, punkt, nelem, elem, nlast, last = lesinput() # -----Plott initalramme----- plot_structure(fig_init, ax_init, punkt, elem, 1, first_index) # -----Regner ut lengder til elementene------ elementlengder = lengder(punkt, elem, nelem) # -----Fastinnspenningsmomentene------ # Lag funksjon selv fim = moment(npunkt, punkt, nelem, elem, nlast, last, elementlengder) # -----Setter opp lastvektor----- # Lag funksjon selv b = lastvektor(fim, npunkt, punkt, nelem, elem) # ------Setter opp systemstivhetsmatrisen----- # Lag funksjon selv K = stivhet(nelem, elem, elementlengder, npunkt) # ------Innfører randbetingelser------ # Lag funksjon selv Kn, Bn = bc(npunkt, punkt, K, b) # -----Løser ligningssystemet------ # Lag funksjon selv rot = ... # Hint, se side for løsing av lineære systemer i Python #------Finner endemoment for hvert element----- # Lag funksjon selv endemoment = endeM(npunkt, punkt, nelem, elem, elementlengder, rot, fim) #-----Skriver ut hva rotasjonen ble i de forskjellige nodene----- print("Rotasjoner i de ulike punktene:") print(rot) #-----Skriver ut hva momentene ble for de forskjellige elementene----- print("Elementvis endemoment:") print(endemoment) #-----Plott deformert ramme----- skalering = 100; # Du kan endre denne konstanten for å skalere de synlige deformasjonene til rammen plot_structure_def(fig_def, ax_def, punkt, elem, 1, first_index, skalering*rot) plt.show() |
Info |
---|
Dersom du ønsker tips til å forbedre kodingen din, besøk siden Tips and tricks for coding. Andre nyttige sider for dette prosjektet kan være LaTeX, Inkscape TexText og Python. |
Utdelte funksjoner
Lengder
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
def lengder(knutepunkt, element, nelem): elementlengder = np.zeros((nelem, 1)) # Beregner elementlengder med Pythagoras' laeresetning for i in range (0, nelem): # OBS! Grunnet indekseringsyntaks i Python-arrays vil ikke denne funksjonen fungere naar vi bare har ett element. dx = knutepunkt[element[i, 0], 0] - knutepunkt[element[i, 1], 0] dy = knutepunkt[element[i, 0], 1] - knutepunkt[element[i, 1], 1] elementlengder[i] = np.sqrt(dx*dx + dy*dy) return elementlengder |
Lesinput
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
def lesinput(): # Åpner inputfilen fid = open("input.txt", "r") # Leser totalt antall punkt npunkt = int(fid.readline()) # 'fid.readline()' leser en linje, 'int(...)' gjør at linjen tolkes som et heltall # LESER INN XY-KOORDINATER TIL KNUTEPUNKTENE # Nodenummer tilsvarer radnummer i "Node-variabel" # x-koordinat lagres i kolonne 1, y-koordinat i kolonne 2 # Grensebetingelse lagres i kolonne 3; 1 = fast innspent og 0 = fri rotasjon punkt = np.loadtxt(fid, dtype = int, max_rows = npunkt) # 'max_rows = npunkt' sorger for at vi bare leser # de 'npunkt' neste linjene i tekstfilen # Leser antall elementer nelem = int(fid.readline()) # Leser konnektivitet: Sammenheng mellom elementender og knutepunktnummer samt EI for elementene # Elementnummer tilsvarer radnummer i "elem"-variabel # Knutepunktnummer for lokal ende 1 lagres i kolonne 1 # Knutepunktnummer for lokal ende 2 lagres i kolonne 2 # Det anbefales at knutepunktnummerering starter på 0, slik at det samsvarerer med listeindeksering i Python # E-modul for materiale lagres i kolonne 3 # Tverrsnittstype lagres i kolonne 4; I-profil = 1 og rørprofil = 2 elem = np.loadtxt(fid, dtype = int, max_rows = nelem) # Leser antall laster som virker på rammen nlast = int(fid.readline()) # Leser lastdata # Bestem selv verdiene som er nødvendig å lese inn, samt hva verdiene som leses inn skal representere last = np.loadtxt(fid, dtype = float, max_rows = nlast) # <-- Forslag til innlesing av lastdata # Lukker input-filen fid.close() return npunkt, punkt, nelem, elem, nlast, last |
Visualisering
Den tildelte visualiseringdelen av prosjektet består av en fil (structure_visualization) med flere funksjoner som er laget for plotte rammen både før og etter lastene er påført. Disse plottene er i hovedsak ikke lagt for å være veldig pene, men det anbefales heller at de blir brukt til å verifisere at input-filen har blitt tolket riktig, samt at de endelige rotasjonene er noenlunde realistiske.
Info |
---|
Det er ikke et krav å bruke disse funksjonene dersom det ikke ønskes, og det er heller ikke et krav å forstå hva som blir gjort inni funksjonene dersom man velger å bruke dem. De er kun laget for å være til hjelp ved å tilby en simpel visualisering av problemet (og delvis løsningen), da dette av erfaring er en mer motiverende måte å debugge på. |
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
# This visualization is based on the original Matlab code by Josef Kiendl, and has been modified to fit TMR4176 by Jon Arnt Kårstad
# NB! Denne filen krever at du har installert Python-pakkene: NumPy, SciPy og Matplotlib
# More detailed information regarding Python (matrices, visualizations etc.) may be found at 'https://www.ntnu.no/wiki/display/imtsoftware'
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline
def setup_plots():
fig_init, ax_init = plt.subplots()
fig_def, ax_def = plt.subplots()
ax_init.set_title('Initialramme')
ax_def.set_title('Deformert ramme')
ax_init.axes.set_aspect('equal')
ax_def.axes.set_aspect('equal')
return fig_init, ax_init, fig_def, ax_def
def plot_structure(fig, ax, punkt, elem, numbers, index_start):
# This is a translation of the original function written by Josef Kiendl in Matlab
# It has been slightly modified in order to be used in TMR4176
# This function plots the beam structure defined by nodes and elements
# The bool (0 or 1) 'numbers' decides if node and element numbers are plotted or not
# Change input to the correct format
nodes = np.array(punkt[:, 0:2], copy = 1, dtype = int)
fixed_dof = np.array(np.where(punkt[:, 2] == 1)[0], copy=1, dtype=int) + 1
el_nod = np.array(elem[:, 0:2], copy=1, dtype=int) + 1
nod_dof = np.arange(1, nodes.shape[0] + 1, 1, dtype=int)
# Start plotting part
for iel in range(0, el_nod.shape[0]):
# Plot element
ax.plot([nodes[el_nod[iel, 0] - 1, 0], nodes[el_nod[iel, 1] - 1, 0]],
[nodes[el_nod[iel, 0] - 1, 1], nodes[el_nod[iel, 1] - 1, 1]], '-k', linewidth = 2)
if numbers == 1:
# Plot element numbers. These are not plotted in the midpoint to
# avoid number superposition when elements cross in the middle
ax.text(nodes[el_nod[iel, 0] - 1, 0] + ( nodes[el_nod[iel, 1] - 1, 0] - nodes[el_nod[iel, 0] - 1, 0] ) / 2.5,
nodes[el_nod[iel, 0] - 1, 1] + ( nodes[el_nod[iel, 1] - 1, 1] - nodes[el_nod[iel, 0] - 1, 1] ) / 2.5,
str(iel + index_start), color = 'blue', fontsize = 16)
if numbers == 1:
# Plot node number
for inod in range(0, nodes.shape[0]):
ax.text(nodes[inod, 0], nodes[inod, 1], str(inod + index_start), color = 'red', fontsize = 16)
def plot_structure_def(fig, ax, punkt, elem, numbers, index_start, r):
# This is a translation of the original function written by Josef Kiendl in Matlab
# This function plots the deformed beam structure defined by nodes and elements
# The bool (0 or 1) 'numbers' decides if node and element numbers are plotted or not
# Change input to the correct format
nodes = np.array(punkt[:, 0:2], copy = 1, dtype = int)
fixed_dof = np.array(np.where(punkt[:, 2] == 1)[0], copy=1, dtype=int) + 1
el_nod = np.array(elem[:, 0:2], copy=1, dtype=int) + 1
nod_dof = np.arange(1, nodes.shape[0] + 1, 1, dtype=int)
if numbers == 1:
# Plot node number
for inod in range(0, nodes.shape[0]):
ax.text(nodes[inod, 0], nodes[inod, 1], str(inod + index_start), color = 'red', fontsize = 16)
for iel in range(0, el_nod.shape[0]):
delta_x = nodes[el_nod[iel, 1] - 1, 0] - nodes[el_nod[iel, 0] - 1, 0]
delta_z = nodes[el_nod[iel, 1] - 1, 1] - nodes[el_nod[iel, 0] - 1, 1]
L = np.sqrt(delta_x ** 2 + delta_z ** 2)
if delta_z >= 0:
psi = np.arccos(delta_x / L)
else:
psi = -np.arccos(delta_x / L)
phi = np.zeros((2, 1))
for inod in range(0, 2):
if nod_dof[el_nod[iel, inod] - 1] > 0:
phi[inod] = r[nod_dof[el_nod[iel, inod] - 1] - 1]
x = np.array([0, L])
z = np.array([0, 0])
xx = np.arange(0, 1.01, 0.01)*L
cs = CubicSpline(x, z, bc_type = ((1, -phi[0, 0]), (1, -phi[1, 0])))
zz = cs(xx)
# Rotate
xxzz = np.array([[np.cos(psi), -np.sin(psi)], [np.sin(psi), np.cos(psi)]]) @ np.vstack([xx, zz])
# Displace
xx2 = xxzz[0, :] + nodes[el_nod[iel, 0] - 1, 0]
zz2 = xxzz[1, :] + nodes[el_nod[iel, 0] - 1, 1]
ax.plot(xx2, zz2, '-k', linewidth = 2)
if numbers == 1:
# Plot element numbers. These are not plotted in the midpoint to
# avoid number superposition when elements cross in the middle
ax.text(xx2[round(xx2.size / 2.5)], zz2[round(xx2.size / 2.5)], str(iel + index_start), color = 'blue', fontsize = 16) |
Inputfil - Eksempel
Følgende er et eksempel på hvordan input.txt kan se ut. I dette eksemplet er det ikke lagt til noen laster på konstruksjonen.
Note |
---|
OBS! Husk å fjerne kommentarene fra egen inputfil. |
No Format |
---|
9 # Antall knutepunkt 0 0 1 # [x, y, Innspenning] 0 10 0 0 20 0 10 20 0 10 10 0 10 0 0 20 20 0 20 10 0 20 0 0 10 # Antall element 0 1 70854000 1 # [Lokal ende 1, Lokal ende 2, Elastisitetsmodul, Profil] 1 2 70854000 1 1 4 141708000 2 2 3 70854000 2 4 3 70854000 1 5 4 70854000 1 3 6 70854000 2 4 7 141708000 1 7 6 70854000 2 8 7 70854000 1 0 # Antall laster |