Moved over from ApproxFun

LICENSE.md

@@ -0,0 +1,22 @@
+The DiskFun.jl package is licensed under the MIT "Expat" License:
+
+> Copyright (c) 2015: Sheehan Olver.
+>
+> Permission is hereby granted, free of charge, to any person obtaining
+> a copy of this software and associated documentation files (the
+> "Software"), to deal in the Software without restriction, including
+> without limitation the rights to use, copy, modify, merge, publish,
+> distribute, sublicense, and/or sell copies of the Software, and to
+> permit persons to whom the Software is furnished to do so, subject to
+> the following conditions:
+>
+> The above copyright notice and this permission notice shall be
+> included in all copies or substantial portions of the Software.
+>
+> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+> EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+> MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+> IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+> CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+> TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+> SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

README.md

@@ -1,2 +1,28 @@
-# DiskFun.jl
-Supports approximating functions and solving differential equations on disks in ApproxFun
+# DiskFun
+
+[![Build Status](https://travis-ci.org/dlfivefifty/DiskFun.jl.svg?branch=master)](https://travis-ci.org/dlfivefifty/DiskFun.jl)
+
+
+
+The following solves Poisson `Δu =f` with zero Dirichlet conditions
+on a disk
+
+```julia
+d = Disk()
+f = Fun((x,y)->exp(-10(x+.2)^2-20(y-.1)^2),d) 
+u = [dirichlet(d);lap(d)]\Any[0.,f]
+ApproxFun.plot(u)                           # Requires Gadfly or PyPlot
+```
+
+
+The following solves beam equation `u_tt + Δ^2u = 0`
+on a disk
+
+```julia
+d = Disk()
+u0 = Fun((x,y)->exp(-50x^2-40(y-.1)^2)+.5exp(-30(x+.5)^2-40(y+.2)^2),d)
+B= [dirichlet(d),neumann(d)]
+L = -lap(d)^2
+h = 0.001
+timeevolution(2,B,L,u0,h)                 # Requires GLPlot
+```

REQUIRE

@@ -0,0 +1 @@
+julia 0.3

appveyor.yml

@@ -0,0 +1,34 @@
+environment:
+  matrix:
+  - JULIAVERSION: "julialang/bin/winnt/x86/0.3/julia-0.3-latest-win32.exe"
+  - JULIAVERSION: "julialang/bin/winnt/x64/0.3/julia-0.3-latest-win64.exe"
+  - JULIAVERSION: "julianightlies/bin/winnt/x86/julia-latest-win32.exe"
+  - JULIAVERSION: "julianightlies/bin/winnt/x64/julia-latest-win64.exe"
+
+branches:
+  only:
+    - master
+    - /release-.*/
+
+notifications:
+  - provider: Email
+    on_build_success: false
+    on_build_failure: false
+    on_build_status_changed: false
+
+install:
+# Download most recent Julia Windows binary
+  - ps: (new-object net.webclient).DownloadFile(
+        $("http://s3.amazonaws.com/"+$env:JULIAVERSION),
+        "C:\projects\julia-binary.exe")
+# Run installer silently, output to C:\projects\julia
+  - C:\projects\julia-binary.exe /S /D=C:\projects\julia
+
+build_script:
+# Need to convert from shallow to complete for Pkg.clone to work
+  - IF EXIST .git\shallow (git fetch --unshallow)
+  - C:\projects\julia\bin\julia -e "versioninfo();
+      Pkg.clone(pwd(), \"DiskFun\"); Pkg.build(\"DiskFun\")"
+
+test_script:
+  - C:\projects\julia\bin\julia --check-bounds=yes -e "Pkg.test(\"DiskFun\")"

src/DiskFun.jl

@@ -0,0 +1,56 @@
+module DiskFun
+    using Base, Compat, ApproxFun
+
+# package code goes here
+import Base: values,getindex,setindex!,*,.*,+,.+,-,.-,==,<,<=,>,
+                >=,./,/,.^,^,\,∪,transpose
+
+importall ApproxFun                
+
+# ApproxFun general import                
+import ApproxFun: BandedMatrix,bazeros,order,
+                  linesum,complexlength,
+                  real, eps, isapproxinteger                
+
+# Operator import    
+import ApproxFun:    bandinds,SpaceOperator, ConversionWrapper, DerivativeWrapper,
+                  rangespace, domainspace, addentries!, BandedOperator, 
+                  promotedomainspace,  CalculusOperator, interlace, Multiplication, 
+                  DiagonalArrayOperator, Recurrence, CompactFunctional, choosedomainspace,
+                    Dirichlet, Neumann, Laplacian, ConstantTimesOperator, Conversion, isfunctional,
+                    dirichlet, neumann, Derivative
+
+
+# Spaces import                  
+import ApproxFun: PolynomialSpace,ConstantSpace, IntervalSpace,
+                    SumSpace,PiecewiseSpace, ArraySpace,RealBasis,ComplexBasis,AnyBasis,
+                    UnsetSpace, AnySpace, canonicalspace, domain, evaluate,
+                    AnyDomain, plan_transform,plan_itransform,
+                    transform,itransform,transform!,itransform!,
+                    isambiguous, fromcanonical, tocanonical, checkpoints, ∂, spacescompatible,
+                   mappoint, UnivariateSpace, setdomain, Space, points, space, conversion_rule, maxspace_rule,
+                   coefficients
+
+# Multivariate import
+import ApproxFun: BivariateDomain,DirectSumSpace,TupleSpace, AbstractProductSpace,    
+                    BivariateFun,  ProductFun, LowRankFun, lap, columnspace, diagop, isproductop, discretize,
+                    schurfact, kronfact, isdiagop
+
+
+# Jacobi import
+import ApproxFun: jacobip, JacobiSD                    
+
+# Plot import
+import ApproxFun: plot, surf
+
+
+
+
+
+
+include("JacobiSquare.jl")
+include("DiskSpace.jl")
+
+include("plot.jl")
+
+end # module

src/DiskSpace.jl

@@ -0,0 +1,192 @@
+
+
+export Disk
+
+##TODO: make argument
+immutable Disk <: BivariateDomain{Float64}
+    radius::Float64
+    center::@compat(Tuple{Float64,Float64})
+end
+
+Disk(r)=Disk(r,(0.,0.))
+Disk()=Disk(1.)
+Disk(::AnyDomain)=Disk(NaN,(NaN,NaN))
+
+
+isambiguous(d::Disk)=isnan(d.radius) && all(isnan,d.center)
+Base.convert(::Type{Disk},::AnyDomain)=Disk(AnyDomain())
+
+
+#canonical is rectangle [r,0]x[-π,π]
+# we assume radius and centre are zero for now
+fromcanonical(D::Disk,x,t)=x*cos(t),x*sin(t)
+tocanonical(D::Disk,x,y)=sqrt(x^2+y^2),atan2(y,x)
+checkpoints(d::Disk)=[fromcanonical(d,(.1,.2243));fromcanonical(d,(-.212423,-.3))]
+
+# function points(d::Disk,n,m,k)
+#     ptsx=0.5*(1-gaussjacobi(n,1.,0.)[1])
+#     ptst=points(PeriodicInterval(),m)
+#
+#     Float64[fromcanonical(d,x,t)[k] for x in ptsx, t in ptst]
+# end
+
+
+∂(d::Disk)=Circle(Complex(d.center...),d.radius)
+
+
+immutable DiskSpace{m,a,b,JS,S} <: AbstractProductSpace{@compat(Tuple{JS,S}),Complex128,2}
+    domain::Disk
+    spacet::S
+    DiskSpace(d,sp)=new(d,sp)
+    DiskSpace(d::AnyDomain)=new(Disk(d),S())
+end
+
+
+DiskSpace(D::Disk,S::FunctionSpace)=DiskSpace{0,0,0,JacobiSquare,typeof(S)}(D,S)
+DiskSpace(D::Disk)=DiskSpace(D,Laurent())
+DiskSpace(d::AnyDomain)=DiskSpace(Disk(d))
+DiskSpace()=DiskSpace(Disk())
+
+canonicalspace(D::DiskSpace)=D
+
+spacescompatible{m,a,b,JS,S}(A::DiskSpace{m,a,b,JS,S},B::DiskSpace{m,a,b,JS,S})=true
+
+coefficient_type{T<:Complex}(::DiskSpace,::Type{T})=T
+coefficient_type{T<:Real}(::DiskSpace,::Type{T})=Complex{T}
+
+domain(d::DiskSpace)=d.domain
+function space(D::DiskSpace,k::Integer)
+    @assert k==2
+    D.spacet
+end
+
+Base.getindex(D::DiskSpace,k::Integer)=space(D,k)
+
+Space(D::Disk)=DiskSpace(D)
+
+
+columnspace{M,a,b,SS}(D::DiskSpace{M,a,b,SS},k)=(m=div(k,2);JacobiSquare(M+m,a+m,b,Interval(D.domain.radius,0.)))
+
+#transform(S::DiskSpace,V::Matrix)=transform([columnspace(S,k) for k=1:size(V,2)],S.spacet,V)
+
+
+evaluate{DS<:DiskSpace}(f::Fun{DS},x...)=evaluate(ProductFun(f),x...)
+
+
+
+function Base.real{JS}(f::ProductFun{JS,Laurent,DiskSpace{0,0,0,JS,Laurent}})
+    cfs=f.coefficients
+    n=length(cfs)
+
+    ret=Array(Fun{JS,Float64},iseven(n)?n+1:n)
+    ret[1]=real(cfs[1])
+
+    for k=2:2:n
+        # exp(1im(k-1)/2*x)=cos((k-1)/2 x) +i sin((k-1)/2 x)
+        ret[k]=imag(cfs[k])
+        ret[k+1]=real(cfs[k])
+    end
+    for k=3:2:n
+        # exp(1im(k-1)/2*x)=cos((k-1)/2 x) +i sin((k-1)/2 x)
+        ret[k]+=real(cfs[k])
+        ret[k-1]-=imag(cfs[k])
+    end
+
+    ProductFun(ret,DiskSpace{0,0,0,JS,Fourier}(space(f).domain,Fourier()))
+end
+#Base.imag{S,T}(u::ProductFun{S,Larent,T})=real(TensorFun(imag(u.coefficients),space(u,2)).').'+imag(TensorFun(real(u.coefficients),space(u,2)).').'
+
+
+
+## Conversion
+# These are placeholders for future
+
+conversion_rule{m,a,b,m2,a2,b2,JS,FS}(A::DiskSpace{m,a,b,JS,FS},
+                                      B::DiskSpace{m2,a2,b2,JS,FS})=DiskSpace{max(m,m2),min(a,a2),min(b,b2),JS,FS}(A.domain,B.spacet)
+
+function coefficients{m,a,b,m2,a2,b2,JS,FS}(cfs::Vector,
+                                            A::DiskSpace{m,a,b,JS,FS},
+                                          B::DiskSpace{m2,a2,b2,JS,FS})
+    g=ProductFun(Fun(cfs,A))
+    rcfs=Fun{typeof(columnspace(B,1)),eltype(cfs)}[Fun(g.coefficients[k],columnspace(B,k)) for k=1:length(g.coefficients)]
+    Fun(ProductFun(rcfs,B)).coefficients
+end
+
+
+# function coefficients{S,V,SS,T}(f::ProductFun{S,V,SS,T},sp::ProductRangeSpace)
+#     @assert space(f,2)==space(sp,2)
+
+#     n=min(size(f,2),length(sp.S))
+#     F=[coefficients(f.coefficients[k],rangespace(sp.S.Rdiags[k])) for k=1:n]
+#     m=mapreduce(length,max,F)
+#     ret=zeros(T,m,n)
+#     for k=1:n
+#         ret[1:length(F[k]),k]=F[k]
+#     end
+#     ret
+# end
+
+
+## Operators
+
+isfunctional{DS<:DiskSpace}(D::Dirichlet{DS},k)=k==1
+
+isdiagop{DS<:DiskSpace}(L::Dirichlet{DS},k)=k==2
+diagop{DS<:DiskSpace}(D::Dirichlet{DS},col)=Evaluation(columnspace(domainspace(D),col),false,D.order)
+
+
+isdiagop{DS<:DiskSpace}(L::Laplacian{DS},k)=k==2
+function diagop{DS<:DiskSpace}(L::Laplacian{DS},col)
+    csp=columnspace(domainspace(L),col)
+    rsp=columnspace(rangespace(L),col)
+    Dt=Derivative(space(domainspace(L),2))
+    c=Dt[col,col]
+
+
+    r=Fun(identity,[domain(L).radius,0.])
+    D=Derivative(csp)
+    Δ=D^2+(1/r)*D+Multiplication((c/r)^2,csp)
+
+    Δ^L.order
+end
+
+
+
+isproductop{DS1<:DiskSpace,DS2<:DiskSpace}(C::Conversion{DS1,DS2})=true
+isdiagop{DS1<:DiskSpace,DS2<:DiskSpace}(C::Conversion{DS1,DS2},k)=k==2
+diagop{DS1<:DiskSpace,DS2<:DiskSpace}(C::Conversion{DS1,DS2},col)=Conversion(columnspace(domainspace(C),col),
+                                                                               columnspace(rangespace(C),col))
+#deprod{DS1<:DiskSpace,DS2<:DiskSpace}(C::Conversion{DS1,DS2},k,::Colon)=ConstantOperator(1.0)
+
+
+lap(d::Disk)=Laplacian(Space(d))
+dirichlet(d::Disk)=Dirichlet(Space(d))
+neumann(d::Disk)=Neumann(Space(d))
+
+lap(d::DiskSpace)=Laplacian(d)
+dirichlet(d::DiskSpace)=Dirichlet(d)
+neumann(d::DiskSpace)=Neumann(d)
+
+
+
+function rangespace{m,a,b,JS,S}(L::Laplacian{DiskSpace{m,a,b,JS,S}})
+    sp=domainspace(L)
+    DiskSpace{m-2L.order,a+2L.order,b+2L.order,JS,S}(sp.domain,sp.spacet)
+end
+
+
+
+# special case of integer modes
+function diagop{b}(L::Laplacian{DiskSpace{0,0,b,JacobiSquare,Laurent}},col)
+    S=columnspace(domainspace(L),col)
+    Dp=DDp(S)
+    Dm=DDm(rangespace(Dp))
+    2Dm*Dp
+end
+
+
+
+function rangespace{b,JS,S}(L::Laplacian{DiskSpace{0,0,b,JS,S}})
+    sp=domainspace(L)
+    DiskSpace{0,0,b+2,JS,S}(sp.domain,sp.spacet)
+end