Abgabe
这个提交包含在:
父节点
d13f1918a0
当前提交
eb4e741ebb
文件差异内容过多而无法显示
加载差异
文件差异内容过多而无法显示
加载差异
文件差异内容过多而无法显示
加载差异
文件差异内容过多而无法显示
加载差异
文件差异内容过多而无法显示
加载差异
文件差异内容过多而无法显示
加载差异
二进制文件未显示。
|
@ -0,0 +1,5 @@
|
|||
Group 1
|
||||
12134031, Bianchi Riccardo
|
||||
01128052, Kapla Daniel
|
||||
01630056, Kuen Jakob
|
||||
01620740, Müller David
|
二进制文件未显示。
二进制文件未显示。
二进制文件未显示。
二进制文件未显示。
二进制文件未显示。
二进制文件未显示。
|
@ -0,0 +1,12 @@
|
|||
module Parameters
|
||||
export k, hz, L, q, T, e, elₘ
|
||||
|
||||
k = 429.
|
||||
hz = 0.0005
|
||||
L = 0.01
|
||||
|
||||
q = 2000000.
|
||||
T = 293.
|
||||
c = 10.
|
||||
elₘ = [41,42,43,44,45,46,47,59,60,61,62,63,77,78,79,95]
|
||||
end
|
|
@ -0,0 +1,6 @@
|
|||
FEM solver
|
||||
|
||||
- programmed in Julia
|
||||
- run using `julia run.jl`
|
||||
- this generates output txt files in ./txt/ and plots in ./plots/ for all 5 variants
|
||||
- the script requires the Julia packages GeometryBasics, GLMakie for plotting, and Fmt (https://github.com/bicycle1885/Fmt.jl) for print_HTP functionality3
|
|
@ -0,0 +1,104 @@
|
|||
module FEM
|
||||
export Node, Element, stiffness, tri_tiles, gradient, center, normalize, p
|
||||
|
||||
include("./Parameters.jl")
|
||||
|
||||
using GeometryBasics
|
||||
using GeometryBasics.LinearAlgebra
|
||||
|
||||
struct Node
|
||||
x::Float64
|
||||
y::Float64
|
||||
index::UInt
|
||||
end
|
||||
|
||||
p(n::Node) = Point2f(n.x, n.y)
|
||||
|
||||
struct Element
|
||||
nodes::Vector{Node}
|
||||
a::Vector{Float64}
|
||||
b::Vector{Float64}
|
||||
c::Vector{Float64}
|
||||
Δ::Float64
|
||||
|
||||
function Element(n::Vector{Node})
|
||||
xs = map(nd->nd.x,n)
|
||||
ys = map(nd->nd.y,n)
|
||||
a = cross(xs, ys)
|
||||
b = cross(ys, ones(3))
|
||||
c = cross(ones(3), xs)
|
||||
Δ = dot(xs,b) / 2
|
||||
|
||||
new(n,a,b,c,Δ)
|
||||
end
|
||||
end
|
||||
|
||||
function stiffness(e::Element)::Matrix{Float64}
|
||||
# tensor product
|
||||
Hₑ = e.b .* e.b'
|
||||
Hₑ += e.c .* e.c'
|
||||
Hₑ *= Parameters.hz*Parameters.k/4e.Δ
|
||||
|
||||
return Hₑ
|
||||
end
|
||||
|
||||
function tri_tiles(L::Float64, divisions::Int, trapezoidal::Bool=false, biased::Bool=false, ring::Bool=false)::Tuple{Vector{Node}, Vector{Element}}
|
||||
nodes = Matrix{Node}(undef, divisions+1, divisions+1)
|
||||
|
||||
i = 1
|
||||
for y in 0:divisions
|
||||
for x in 0:divisions
|
||||
xₑ = L*x/divisions
|
||||
yₑ = L*y/divisions
|
||||
|
||||
if biased
|
||||
B = yₑ/2Parameters.L
|
||||
xₑ = xₑ*(xₑ*B/Parameters.L - B + 1)
|
||||
end
|
||||
|
||||
if trapezoidal
|
||||
xₑ *= 1 - 0.5*(y/divisions)
|
||||
end
|
||||
|
||||
if ring
|
||||
r = Parameters.L*(1 + x/divisions)
|
||||
θ = (π/4)*(y/divisions)
|
||||
|
||||
xₑ = r*cos(θ)
|
||||
yₑ = r*sin(θ)
|
||||
xₑ -= Parameters.L
|
||||
end
|
||||
|
||||
nodes[x+1,y+1] = Node(xₑ,yₑ,i)
|
||||
i += 1
|
||||
end
|
||||
end
|
||||
|
||||
elements = []
|
||||
|
||||
for y in 1:divisions
|
||||
for x in 1:divisions
|
||||
# lower/upper triangle
|
||||
push!(elements, Element([nodes[x,y], nodes[x+1,y], nodes[x,y+1]]))
|
||||
push!(elements, Element([nodes[x+1,y+1], nodes[x,y+1], nodes[x+1,y]]))
|
||||
end
|
||||
end
|
||||
|
||||
return vec(nodes), elements
|
||||
end
|
||||
|
||||
function gradient(e::Element, T::Vector{Float64})::Vec2f
|
||||
return Vec2f([
|
||||
e.b[1] e.b[2] e.b[3] ;
|
||||
e.c[1] e.c[2] e.c[3]
|
||||
] * [
|
||||
T[e.nodes[1].index]
|
||||
T[e.nodes[2].index]
|
||||
T[e.nodes[3].index]
|
||||
] / 2e.Δ)
|
||||
end
|
||||
|
||||
center(e::Element)::Point2f = Point2f(sum([n.x for n in e.nodes])/3.0, sum([n.y for n in e.nodes])/3.0)
|
||||
normalize(v::Vec2f)::Vec2f = v / sqrt(v[1]^2 + v[2]^2)
|
||||
|
||||
end
|
|
@ -0,0 +1,39 @@
|
|||
using Fmt
|
||||
|
||||
function print_HTP(H::Matrix{Float64}, T::Vector{Float64}, P::Vector{Float64}, filename="output.txt")
|
||||
# Print matrices to .txt-file (name of file = filename).
|
||||
# H... overall assembled stiffness matrix
|
||||
# T... nodal temperature vector
|
||||
# P... nodal force vector
|
||||
|
||||
# Make sure, that your system of equations is sorted by
|
||||
# ascending node numbers, i.e., N1 N2 ... N100.
|
||||
|
||||
open(filename, "w") do io
|
||||
write(io, "Stiffness matrix H: \n")
|
||||
|
||||
for row in H
|
||||
for col in row
|
||||
outline = f"{$col:+8.4e},"
|
||||
write(io, f"{$outline:11s}")
|
||||
end
|
||||
write(io, "\n")
|
||||
end
|
||||
|
||||
write(io, "Temperature T: \n")
|
||||
for row in T
|
||||
for col in row
|
||||
outline = f"{$col:+8.4e},"
|
||||
write(io, f"{$outline:11s} \n")
|
||||
end
|
||||
end
|
||||
|
||||
write(io, "Force vector P: \n")
|
||||
for row in P
|
||||
for col in row
|
||||
outline = f"{$col:+8.4e},"
|
||||
write(io, f"{$outline:11s} \n")
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
|
@ -0,0 +1,94 @@
|
|||
include("./fem.jl")
|
||||
include("./Parameters.jl")
|
||||
include("./print_HTP.jl")
|
||||
|
||||
using .FEM
|
||||
using GLMakie
|
||||
using Makie.GeometryBasics
|
||||
|
||||
@enum Variation BASIC V1 V2 V3 V4_div V4_mul
|
||||
|
||||
function solve(name::String, variation::Variation)
|
||||
div = 9
|
||||
H = zeros(Float64,(div+1)^2, (div+1)^2)
|
||||
|
||||
(nodes, elements) = tri_tiles(Parameters.L, div, variation==V1, variation==V2, variation==V3)
|
||||
|
||||
for (index, el) in enumerate(elements)
|
||||
Hₑ = stiffness(el)
|
||||
|
||||
if variation == V4_div && index in Parameters.elₘ
|
||||
Hₑ /= Parameters.c
|
||||
elseif variation == V4_mul && index in Parameters.elₘ
|
||||
Hₑ *= Parameters.c
|
||||
end
|
||||
|
||||
for i in 1:3
|
||||
for j in 1:3
|
||||
H[el.nodes[i].index, el.nodes[j].index] += Hₑ[i,j]
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
P = zeros(Float64, 90)
|
||||
|
||||
# impose neumann conditions on top edge
|
||||
# adapted from https://mathoverflow.net/questions/5085/how-to-apply-neuman-boundary-condition-to-finite-element-method-problems
|
||||
# maybe check if this actually makes sense
|
||||
|
||||
nᵧ = -Parameters.q / Parameters.k
|
||||
# skip=2 since only every second element has a edge along the bottom: ◺◹
|
||||
for ∂_el in elements[1:2:18]
|
||||
n1 = ∂_el.nodes[1]
|
||||
n2 = ∂_el.nodes[2]
|
||||
l = abs(n2.x - n1.x)
|
||||
P[n1.index] += nᵧ*l/2
|
||||
P[n2.index] += nᵧ*l/2
|
||||
end
|
||||
|
||||
rhs = P - H[1:90,91:100]*fill(Parameters.T, 10,1)
|
||||
T = vec(H[1:90,1:90]\rhs)
|
||||
append!(T, fill(Parameters.T, 10))
|
||||
|
||||
reaction_forces = H[91:100,:]*T
|
||||
|
||||
centers = center.(elements)
|
||||
|
||||
gradients = map(el -> gradient(el, T), elements)
|
||||
flux = gradients .* -1
|
||||
|
||||
norm = maximum(map(g-> sqrt(g[1]^2 + g[2]^2), gradients))/(Parameters.L/div)*2.5
|
||||
|
||||
gradients ./= norm
|
||||
flux ./= norm
|
||||
|
||||
set_theme!(theme_black())
|
||||
f = Figure(resolution = (1536, 1024))
|
||||
|
||||
tris = map(e->Polygon([p(e.nodes[1]), p(e.nodes[2]), p(e.nodes[3])]), elements)
|
||||
poly(f[1, 1], tris, color=:transparent, linestyle=:solid, strokewidth=0.8, strokecolor=:white, transparency=true)
|
||||
|
||||
xs = map(n->n.x, nodes)
|
||||
ys = map(n->n.y, nodes)
|
||||
|
||||
xs = reshape(xs, (div+1, div+1))
|
||||
ys = reshape(ys, (div+1, div+1))
|
||||
Ts = reshape(T, (div+1, div+1))
|
||||
|
||||
surface(f[2, 1], xs, ys, Ts, colormap=:matter, axis=(type=Axis3,))
|
||||
|
||||
contour(f[1:2,2:3], map(n->n.x,nodes), map(n->n.y,nodes), T, levels=16, colormap=:matter)
|
||||
arrows!(f[1:2,2:3], centers, gradients, arrowcolor=:red, linecolor=:red)
|
||||
arrows!(f[1:2,2:3], centers, flux, arrowcolor=:blue, linecolor=:blue)
|
||||
|
||||
save("plots/$name.png", current_figure())
|
||||
|
||||
# print_HTP(H, T, P, "txt/htp_$name.txt")
|
||||
end
|
||||
|
||||
solve("baseline", BASIC)
|
||||
solve("variation1", V1)
|
||||
solve("variation2", V2)
|
||||
solve("variation3", V3)
|
||||
solve("variation4_div", V4_div)
|
||||
solve("variation4_mul", V4_mul)
|
正在加载...
在新工单中引用