Skip to content

Commit

Permalink
upgraded to TensorField methods
Browse files Browse the repository at this point in the history
  • Loading branch information
@chakravala
chakravala committed Nov 27, 2024
1 parent 3e03293 commit 07f7e22
Show file tree
Hide file tree
Showing 4 changed files with 168 additions and 81 deletions.
4 changes: 2 additions & 2 deletions Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "FlowGeometry"
uuid = "af326652-bfa1-11e9-1b34-cf9047ce1632"
authors = ["Michael Reed"]
version = "0.1.0"
version = "0.1.1"

[deps]
Cartan = "e24af43f-9fd5-4672-9264-a75f3ae40eb2"
Expand All @@ -12,7 +12,7 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df"
[compat]
julia = "1"
Cartan = "0.3"
Grassmann = "0.8, 0.9"
Grassmann = "0.8"
Requires = "1"

[extras]
Expand Down
196 changes: 134 additions & 62 deletions src/FlowGeometry.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,8 @@ function edgeslist!(p,r)
Values{2,Int}.(l+1:l+n,l.+[2:n;1])
end

# need PointCloud etc operations for refinement (same ID) or union etc (new ID)

addbound(e,r=rectangle(xn,xm,yn,ym)) = [e;edgeslist!(points(e),r)]
airfoilbox(n) = addbound(airfoiledges(n),rectangle(xn,xm,yn,ym))
airfoiledges(n) = edgeslist(PointCloud(points(n)))
Expand Down Expand Up @@ -81,15 +83,20 @@ end

export rectangletriangle, rectangletriangles, FittedPoint

rectangletriangle(i,m,p) = rectangletriangle(((i-1)%(2(m-1)))+1,((i-1)÷(2(m-1)))+1,m,p)
function rectangletriangle(i,j,m,p)
rectangletriangle(i,m) = rectangletriangle(((i-1)%(2(m-1)))+1,((i-1)÷(2(m-1)))+1,m)
function rectangletriangle(i,j,m)
k=(j-1)*m+(i÷2)+1;n=m+k
isodd(i) ? Values(k,n,k+1) : Values(k,n-1,n)
end
#rectangletriangle(i,j,m,p) = (k=(j-1)*m+(i÷2)+1;n=m+k;Chain{p,1}(iseven(i) ? Values(k,n,n+1) : Values(k,k-1,n)))
#rectangletriangle(i,j,m) = (k=(j-1)*m+(i÷2)+1;n=m+k; iseven(i) ? Values(k,n,n+1) : Values(k,k-1,n))

rectangletriangles(p,m=51,JL=51) = [rectangletriangle(i,JL,p) for i ∈ 1:2*(m-1)*(JL-1)]
rectanglebounds(n=51,JL=51) = [1:JL:JL*n; JL*(n-1)+2:JL*n; JL*(n-1):-JL:JL; JL-1:-1:2]
function rectangletriangles(m=51,JL=51)
SimplexManifold([rectangletriangle(i,JL) for i ∈ 1:2*(m-1)*(JL-1)],m*JL)
end
function rectanglebounds(n=51,JL=51)
b = [1:JL:JL*n; JL*(n-1)+2:JL*n; JL*(n-1):-JL:JL; JL-1:-1:2]; push!(b,1)
SimplexManifold(Values{2,Int}.([(b[i],b[i+1]) for i ∈ 1:length(b)-1]),n*JL)
end

FittedPoint(k,JL=51) = Chain{Submanifold(ℝ^3),1}(1.0,(k-1)÷JL,(k-1)%JL)

Expand Down Expand Up @@ -120,20 +127,17 @@ function RakichPoint(k,x,y,s,κ,JL=legnth(y))
Chain{Submanifold(ℝ^3),1}(1.0,x[xk],s[xk]≠0 ? Rakich(κ[xk],yk,s[xk],y[end],JL) : y[yk])
end

function RakichPoints(P::CircularArc{T,m}=CircularArc{6,21}(),D=50,n=51,JL=51) where {T,m}
function rakichpoints(P::CircularArc{T,m}=CircularArc{6,21}(),D=50,n=51,JL=51) where {T,m}
t,c,x0 = T/100,1,0
x = RakichPlate(P,D,n)
y = RakichLine(0,D,JL,t/10)
s = profile.(Ref(P),x,c,x0)
κ = [k≠0 ? RakichNewton(D-k,JL,t/10) : 0.0 for k ∈ s]
[RakichPoint(k,x,y,s,κ,JL) for k ∈ 1:n*JL]
PointCloud([RakichPoint(k,x,y,s,κ,JL) for k ∈ 1:n*JL])
end

function initrakich(P=CircularArc{6,61}(),D=50,n=101,JL=51)
p = PointCloud(RakichPoints(P,D,n,JL))
b = rectanglebounds(n,JL); push!(b,1)
e = SimplexManifold(Values{2,Int}.([(b[i],b[i+1]) for i ∈ 1:length(b)-1]),length(p))
return p,e,SimplexManifold(rectangletriangles(p,n),length(p))
rakichpoints(P,D,n,JL),rectanglebounds(n,JL),rectangletriangles(n,JL)
end

# points
Expand Down Expand Up @@ -186,78 +190,146 @@ function sphere(fac::Vector,r=1,p=points(fac))
return out
end

function wing(N::Airfoil,λ=0.7,σ=0.5)
U,L = imag.(upper(N)),imag.(lower(N))
np = length(U)
out1 = zeros(np,2np-1)
out2 = zeros(np,2np-1)
out3 = zeros(np,2np-1)
int = interval(N)
rint = reverse(int)
pf = profile(N.c)
for i ∈ 1:np-1
out1[:,i] .= (σ*int[i]).+(λ+(1-λ)*rint[i]).*int
out2[:,i] .= int[i]
out3[:,i] .= fiber(rint[i]*U)
end
out1[:,np] .= (σ*int[np]).+(λ*rint[np]).*int
out2[:,np] .= int[np]
for i ∈ 1:np-1
out1[:,i+np] .= (σ*rint[i+1]).+(λ+(1-λ)*int[i+1]).*int
out2[:,i+np] .= rint[i+1]
out3[:,i+np] .= fiber(int[i+1]*L)
end
TensorField(int⊕(-1:step(int):1),Chain.(out1,out2,out3))
end

function convhull(p::PointCloud)
n = length(p)
out = Values{2,Int}[]
for i ∈ 1:n
pi = base(p)[i]
for j ∈ 1:n
i==j && continue
pj = base(p)[j]
pij = pi∧pj
val = true
for k ∈ 1:n
k∈(i,j) && continue
pk = base(p)[k]
pijk = pij∧pk
if Real(pijk) ≈ 0
if Real(abs((pi+pj)/2-pk)) < Real(abs(pi-pj))
val = false
break
end
elseif Real(pijk) > 0
val = false
break
end
end
val && push!(out,Values(i,j))
end
end
return SimplexManifold(out,n)
end

function convhull(p::PointCloud,r)
n = length(p)
out = Values{2,Int}[]
for i ∈ 1:n
pi = base(p)[i]
for j ∈ 1:n
i==j && continue
pj = base(p)[j]
pij = pi∧pj
avg = (pi+pj)/2
val = true
for k ∈ 1:n
k∈(i,j) && continue
pk = base(p)[k]
dist = Real(abs(avg-pk))
dist > r && continue
pijk = pij∧pk
if Real(pijk) ≈ 0
if dist < Real(abs(pi-pj))
val = false
break
end
elseif Real(pijk) > 0
val = false
break
end
end
val && push!(out,Values(i,j))
end
end
return SimplexManifold(out,n)
end


# init

function __init__()
@require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" begin
function Makie.plot(N::Profile,args...)
Makie.plot(interval(N),profile(N),args...)
end
function Makie.plot!(N::Profile,args...)
Makie.plot!(interval(N),profile(N),args...)
end
function Makie.plot(N::Airfoil,args...)
Makie.lines(N,args...)
Makie.plot!(N.c,args...)
end
function Makie.plot!(N::Airfoil,args...)
Makie.lines!(N,args...)
Makie.plot!(N.c,args...)
end
Makie.lines(N::Profile,args...) = Makie.lines(profile(N),args...)
Makie.lines!(N::Profile,args...) = Makie.lines!(profile(N),args...)
function Makie.lines(N::Airfoil,args...)
P = complex(N)
Makie.lines([real.(P) imag.(P)],args...)
display(Makie.lines(N,args...))
Makie.lines!(N.c,args...)
end
function Makie.lines!(N::Airfoil,args...)
P = complex(N)
Makie.lines!([real.(P) imag.(P)],args...)
Makie.lines!(N,args...)
Makie.lines!(N.c,args...)
end
function Makie.plot(N::DoubleArc,args...)
Makie.lines(N::Airfoil,args...) = Makie.lines(complex(N),args...)
Makie.lines!(N::Airfoil,args...) = Makie.lines!(complex(N),args...)
function Makie.lines(N::DoubleArc,args...)
U,L = upperlower(N)
Makie.plot(real.(U),imag.(U),args...)
length(U)==length(L) && Makie.plot!(real.(U),(imag.(U).+imag.(L))./2,args...)
Makie.plot!(real.(L),imag.(L),args...)
display(Makie.lines(U,args...))
length(U)==length(L) && Makie.lines!(fiber(real.(U)),fiber(imag.(U).+imag.(L))./2,args...)
Makie.lines!(L,args...)
end
function Makie.plot!(N::DoubleArc,args...)
function Makie.lines!(N::DoubleArc,args...)
U,L = upperlower(N)
Makie.plot!(real.(U),imag.(U),args...)
length(U)==length(L) && Makie.plot!(real.(U),(imag.(U).+imag.(L))./2,args...)
Makie.plot!(real.(L),imag.(L),args...)
Makie.lines!(U,args...)
length(U)==length(L) && Makie.lines!(fiber(real.(U)),fiber(imag.(U).+imag.(L))./2,args...)
Makie.lines!(L,args...)
end
end
@require UnicodePlots="b8865327-cd53-5732-bb35-84acbb429228" begin
function UnicodePlots.Plot(N::Profile,args...)
UnicodePlots.lineplot(interval(N),profile(N),args...)
end
function Plot!(p,N::Profile,args...)
UnicodePlots.lineplot!(p,interval(N),profile(N),args...)
end
function UnicodePlots.Plot(N::Airfoil,args...)
Plot!(UnicodePlots.lineplot(N,args...),N.c,args...)
end
function Plot!(p,N::Airfoil,args...)
Plot!(UnicodePlots.lineplot!(p,N,args...),N.c,args...)
end
UnicodePlots.lineplot(N::Profile,args...) = UnicodePlots.lineplot(profile(N),args...)
UnicodePlots.lineplot!(p,N::Profile,args...) = UnicodePlots.lineplot!(p,profile(N),args...)
function UnicodePlots.lineplot(N::Airfoil,args...)
P = complex(N)
UnicodePlots.lineplot(real.(P),imag.(P),args...)
UnicodePlots.lineplot!(UnicodePlots.lineplot(complex(N),args...),N.c,args...)
end
function UnicodePlots.lineplot!(p,N::Airfoil,args...)
P = complex(N)
UnicodePlots.lineplot!(p,real.(P),imag.(P),args...)
UnicodePlots.lineplot!(UnicodePlots.lineplot!(p,complex(N),args...),N.c,args...)
end
function UnicodePlots.Plot(N::DoubleArc,args...)
function UnicodePlots.lineplot(N::DoubleArc,args...)
U,L = upperlower(N)
p = UnicodePlots.lineplot(real.(U),imag.(U),args...)
length(U)==length(L) && UnicodePlots.lineplot!(p,real.(U),(imag.(U).+imag.(L))./2,args...)
UnicodePlots.lineplot!(p,real.(L),imag.(L),args...)
p = UnicodePlots.lineplot(U,args...)
length(U)==length(L) && UnicodePlots.lineplot!(p,fiber(real.(U)),fiber(imag.(U).+imag.(L))./2,args...)
UnicodePlots.lineplot!(p,L,args...)
end
function Plot!(p,N::DoubleArc,args...)
function UnicodePlots.lineplot!(p,N::DoubleArc,args...)
U,L = upperlower(N)
UnicodePlots.lineplot!(p,real.(U),imag.(U),args...)
length(U)==length(L) && UnicodePlots.lineplot!(p,real.(U),(imag.(U).+imag.(L))./2,args...)
UnicodePlots.lineplot!(p,real.(L),imag.(L),args...)
UnicodePlots.lineplot!(p,U,args...)
length(U)==length(L) && UnicodePlots.lineplot!(p,fiber(real.(U)),fiber(imag.(U).+imag.(L))./2,args...)
UnicodePlots.lineplot!(p,L,args...)
end
Base.display(t::Airfoil) = (display(typeof(t)); display(UnicodePlots.lineplot(t)))
Base.display(t::Profile) = (display(typeof(t)); display(UnicodePlots.lineplot(t)))
end
@require TetGen="c5d3f3f7-f850-59f6-8a2e-ffc6dc1317ea" begin
spheremesh() = TetGen.tetrahedralize(cubesphere(),"vpq1.414a0.1";holes=[Point(0.0,0.0,0.0)])
Expand All @@ -274,7 +346,7 @@ function __init__()
@require MATLAB="10e44e05-a98a-55b3-a45b-ba969058deb6" begin
decsg(args...) = MATLAB.mxcall(:decsg,1,args...)
function decsg(N::Airfoil{pts}) where pts
P = complex(N)[1:end-1]
P = fiber(complex(N))[1:end-1]
x1,x2,x3,x4 = -1.5,3.5,3.5,-1.5
y1,y2,y3,y4 = 1.5,1.5,-1.5,-1.5
R = [[3,4,x1,x2,x3,x4,y1,y2,y3,y4];zeros(2pts-8)]
Expand Down
34 changes: 21 additions & 13 deletions src/airfoils.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,31 +34,31 @@ for Side ∈ ("Upper","Lower")
@pure $Arc{A,p}(a) where {A,p} = new{A,p}(a)
end
@pure $Arc(s::S) where S<:Airfoil{P} where P = $Arc{S,length($side(s))}(s)
$side(z::A,c=1,x0=0) where A<:Profile = $side(interval(z,c,x0),profile(z).*(c*im))
$side(z::A,c=1,x0=0) where A<:Profile = $side(interval(z,c,x0),profile(z)*(c*im))
interval(U::$Arc,c=1,x0=0) = real.($side(U.a,c,x0))
profile(U::$Arc) = imag.($side(U.a))
profileslope(f::$Arc,e=initedges(f)) = interp(e,column(gradient(e,profile(f))))
profileslope(f::$Arc) = gradient(profile(f))
end
end

@pure chord(p) = range(0,0,length=p)
@pure chord(::Profile{p}) where p = chord(p)
@pure chord(::Airfoil{p}) where p = chord(p)
upper(z,r) = (U = z.+r; U[end] = 1; U) # upper surface points xu+im*yu
lower(z,r) = (L = z.-r; L[end] = 1; L) # lower surface points xl+im*yl
upper(x,yc,dyc_dx,yt) = upper(x.+im.*yc,im.*cis.(atan.(dyc_dx)).*yt)
lower(x,yc,dyc_dx,yt) = lower(x.+im.*yc,im.*cis.(atan.(dyc_dx)).*yt)
upperlower(x,yc,dyc_dx,yt) = upperlower(x.+im.*yc,im.*cis.(atan.(dyc_dx)).*yt)
upperlower(z::A,c=1,x=0) where A<:Profile = upperlower(interval(z,c,x),profile(z).*(c*im))
upper(z,r) = (U = z+r; U[end] = 1; U) # upper surface points xu+im*yu
lower(z,r) = (L = z-r; L[end] = 1; L) # lower surface points xl+im*yl
upper(x,yc,dyc_dx,yt) = upper(x+im*yc,im*cis.(atan.(dyc_dx))*yt)
lower(x,yc,dyc_dx,yt) = lower(x+im*yc,im*cis.(atan.(dyc_dx))*yt)
upperlower(x,yc,dyc_dx,yt) = upperlower(x+im*yc,im*cis.(atan.(dyc_dx))*yt)
upperlower(z::A,c=1,x=0) where A<:Profile = upperlower(interval(z,c,x),profile(z)*(c*im))
upperlower(z::A,c=1,x=0) where A<:Airfoil = upper(z,c,x),lower(z,c,x)
upperlower(z,r) = upper(z,r),lower(z,r)

function Base.complex(N::A) where A<:Airfoil
U,L = upperlower(N)
[U;reverse(L)[2:end]]
TorusTopology(TensorField(doubleinterval(interval(N)),[fiber(U);reverse(fiber(L))[2:end]]))
end
function points(N::A) where A<:Airfoil
z = complex(N)[1:end-1]
z = fiber(complex(N))[1:end-1]
Chain{Submanifold(ℝ^3),1}.(1.0,real.(z),imag.(z))
end

Expand Down Expand Up @@ -114,6 +114,8 @@ end
@pure British{n,p}() where {n,p} = British(Camber{Int(n÷100),p}(),ClarkY{n%100,p}())
@pure British(c::C,t::T) where {C<:Profile{P},T<:Profile{P}} where P = British{C,T,2P-2}(c,t)

camber(n::British) = n.c
thickness(n::British) = n.t
upper(n::British,c=1,x0=0) = upper(getprofiles(n,c,x0)...)
lower(n::British,c=1,x0=0) = lower(getprofiles(n,c,x0)...)
upperlower(n::British,c=1,x0=0) = upperlower(getprofiles(n,c,x0)...)
Expand All @@ -135,10 +137,13 @@ end
@pure American{n,p}() where {n,p} = American(Camber{Int(n÷100),p}(),ClarkY{n%100,p}())
@pure American(c::C,t::T) where {C<:Profile{P},T<:Profile{P}} where P = American{C,T,2P-2}(c,t)

camber(n::American) = n.c
thickness(n::American) = n.t
upper(n::American,c=1,x0=0) = upper(getprofiles(n,c,x0)...)
lower(n::American,c=1,x0=0) = lower(getprofiles(n,c,x0)...)
upperlower(n::American,c=1,x0=0) = upperlower(getprofiles(n,c,x0)...)
getprofiles(n::American,c=1,x0=0) = interval(n.c,c,x0),c*profile(n.c),profileslope(n.c),c*profile(n.t)
getprofiles(n::American,c=1,x0=0) = interval(n,c,x0),c*profile(n.c),profileslope(n.c),c*profile(n.t)
interval(n::American,c=1,x0=0) = interval(n.c,c,x0)

# NACA specification selection

Expand Down Expand Up @@ -169,7 +174,10 @@ end

joukowski(R,f,g,b,p) = Joukowski{R,f,g,b,p}()

interval(::Joukowski{R,f,g,b,p}) where {R,f,g,b,p} = interval(2p-1,2π)

function Base.complex(::Joukowski{R,f,g,b,p}) where {R,f,g,b,p}
z = R*cis.(interval(2p-1,2π)).-(f-g*im)
z .+ b^2*inv.(z)
i = interval(2p-1,2π)
z = R*cis.(i).-(f-g*im)
TensorField(i, z .+ b^2*inv.(z))
end
15 changes: 11 additions & 4 deletions src/profiles.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,17 +25,24 @@ export NACA4, NACA5, NACA6, NACA6A

abstract type Profile{P} end

profile(p::T) where T<:Profile = p.(interval(p))
profileslope(p::T) where T<:Profile = profileslope.(Ref(p),interval(p))
function profile(p::T) where T<:Profile
int = interval(p)
TensorField(int,p.(int))
end
function profileslope(p::T) where T<:Profile
int = interval(p)
TensorField(int,profileslope.(Ref(p),int))
end
profileangle(p::T) where T<:Profile = atan.(profileslope(p))
@pure profile(p::T,x) where T<:Profile = p(x)
@pure profile(p::T,x,c,x0=0) where T<:Profile = profile(p,(x-x0)/c)*c
@pure profileslope(p::T,x,c,x0=0) where T<:Profile = profileslope(p,(x-x0)/c)
@pure profileangle(p::T,x,c,x0=0) where T<:Profile = atan(profileslope(p,x,c,x0))
@pure interval(::P,c=1,x0=0) where P<:Profile{p} where p = interval(p,c,x0)
@pure interval(p::Int,c=1,x0=0) = range(x0,c,length=p)
doubleinterval(p) = p[1]:step(p):p[1]+2(p[end]-p[1])
function points(N::A,c=1,x0=0) where A<:Profile{p} where p
Chain{Submanifold(ℝ^3),1}.(1.0,interval(N,c,x0),profile(N))
Chain{Submanifold(ℝ^3),1}.(1.0,interval(N,c,x0),fiber(profile(N)))
end

const AppendixI = [0,.005,.0125,.025,.05,.075,.1,.15,.2,.25,.3,.4,.5,.6,.7,.8,.9,.95,1]
Expand All @@ -55,7 +62,7 @@ end
@pure (n::FlatPlate)(x,c,x0=0) = profile(n,x,c,x0)
@pure profile(::FlatPlate,x,c,x0=nothing) = 0.0
@pure profileslope(::FlatPlate,x,c=nothing,x0=nothing) = 0.0
@pure profile(f::FlatPlate{p}) where p = range(0,0,length=p)
@pure profile(f::FlatPlate{p}) where p = TensorField(interval(f),range(0,0,length=p))
@pure profileslope(f::FlatPlate) = profile(f)

# parabolic arc
Expand Down

2 comments on commit 07f7e22

@chakravala
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register()

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/120307

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.1.1 -m "<description of version>" 07f7e22e15492bd3d4b04d8314312e0d65708a9b
git push origin v0.1.1

Please sign in to comment.