|
1 | 1 | export ContinuousSpace |
2 | 2 |
|
3 | | -struct ContinuousSpace{T,R,P<:PiecewiseSegment{T}} <: Space{P,R} |
4 | | - domain::P |
5 | | -end |
6 | | - |
7 | | -ContinuousSpace(d::PiecewiseSegment{T}) where {T} = |
8 | | - ContinuousSpace{T,real(eltype(T)),typeof(d)}(d) |
9 | | - |
10 | | - |
11 | 3 | Space(d::PiecewiseSegment) = ContinuousSpace(d) |
12 | 4 |
|
13 | | -isperiodic(C::ContinuousSpace) = isperiodic(domain(C)) |
14 | | - |
15 | 5 | const PiecewiseSpaceReal{CD} = PiecewiseSpace{CD,<:Any,<:Real} |
16 | 6 | const PiecewiseSpaceRealChebyDirichlet11 = |
17 | 7 | PiecewiseSpaceReal{<:TupleOrVector{ChebyshevDirichlet{1,1}}} |
18 | 8 |
|
19 | | -spacescompatible(a::ContinuousSpace,b::ContinuousSpace) = domainscompatible(a,b) |
20 | 9 | conversion_rule(a::ContinuousSpace, b::PiecewiseSpaceRealChebyDirichlet11) = a |
21 | 10 |
|
22 | | -plan_transform(sp::ContinuousSpace,vals::AbstractVector) = |
23 | | - TransformPlan{eltype(vals),typeof(sp),false,Nothing}(sp,nothing) |
24 | | - |
25 | | -function *(P::TransformPlan{T,<:ContinuousSpace,false},vals::AbstractVector{T}) where {T} |
26 | | - S = P.space |
27 | | - n=length(vals) |
28 | | - d=domain(S) |
29 | | - K=ncomponents(d) |
30 | | - k=div(n,K) |
31 | | - |
32 | | - PT=promote_type(prectype(d),eltype(vals)) |
33 | | - if k==0 |
34 | | - vals |
35 | | - elseif isperiodic(d) |
36 | | - ret=Array{PT}(undef, max(K,n-K)) |
37 | | - r=n-K*k |
38 | | - |
39 | | - for j=1:r |
40 | | - cfs=transform(ChebyshevDirichlet{1,1}(component(d,j)), |
41 | | - vals[(j-1)*(k+1)+1:j*(k+1)]) |
42 | | - if j==1 |
43 | | - ret[1]=cfs[1]-cfs[2] |
44 | | - ret[2]=cfs[1]+cfs[2] |
45 | | - elseif j < K |
46 | | - ret[j+1]=cfs[1]+cfs[2] |
47 | | - end |
48 | | - ret[K+j:K:end]=cfs[3:end] |
49 | | - end |
50 | | - |
51 | | - for j=r+1:K |
52 | | - cfs=transform(ChebyshevDirichlet{1,1}(component(d,j)), |
53 | | - vals[r*(k+1)+(j-r-1)*k+1:r*(k+1)+(j-r)*k]) |
54 | | - if length(cfs)==1 && j <K |
55 | | - ret[j+1]=cfs[1] |
56 | | - elseif j==1 |
57 | | - ret[1]=cfs[1]-cfs[2] |
58 | | - ret[2]=cfs[1]+cfs[2] |
59 | | - elseif j < K |
60 | | - ret[j+1]=cfs[1]+cfs[2] |
61 | | - end |
62 | | - ret[K+j:K:end]=cfs[3:end] |
63 | | - end |
64 | | - |
65 | | - ret |
66 | | - else |
67 | | - ret=Array{PT}(undef, n-K+1) |
68 | | - r=n-K*k |
69 | | - |
70 | | - for j=1:r |
71 | | - cfs=transform(ChebyshevDirichlet{1,1}(component(d,j)), |
72 | | - vals[(j-1)*(k+1)+1:j*(k+1)]) |
73 | | - if j==1 |
74 | | - ret[1]=cfs[1]-cfs[2] |
75 | | - end |
76 | | - |
77 | | - ret[j+1]=cfs[1]+cfs[2] |
78 | | - ret[K+j+1:K:end]=cfs[3:end] |
79 | | - end |
80 | | - |
81 | | - for j=r+1:K |
82 | | - cfs=transform(ChebyshevDirichlet{1,1}(component(d,j)), |
83 | | - vals[r*(k+1)+(j-r-1)*k+1:r*(k+1)+(j-r)*k]) |
84 | | - |
85 | | - if length(cfs) ≤ 1 |
86 | | - ret .= cfs |
87 | | - else |
88 | | - if j==1 |
89 | | - ret[1]=cfs[1]-cfs[2] |
90 | | - end |
91 | | - ret[j+1]=cfs[1]+cfs[2] |
92 | | - ret[K+j+1:K:end]=cfs[3:end] |
93 | | - end |
94 | | - end |
95 | | - |
96 | | - ret |
97 | | - end |
98 | | -end |
99 | | - |
100 | 11 | components(S::ContinuousSpace) = [ChebyshevDirichlet{1,1}(s) for s in components(domain(S))] |
101 | | -canonicalspace(S::ContinuousSpace) = PiecewiseSpace(components(S)) |
102 | | -convert(::Type{PiecewiseSpace}, S::ContinuousSpace) = canonicalspace(S) |
103 | | - |
104 | | -blocklengths(C::ContinuousSpace) = Fill(ncomponents(C.domain),∞) |
105 | | - |
106 | | -block(C::ContinuousSpace,k) = Block((k-1)÷ncomponents(C.domain)+1) |
107 | | - |
108 | | - |
109 | | -## components |
110 | | - |
111 | | -components(f::Fun{<:ContinuousSpace},j::Integer) = components(Fun(f,canonicalspace(f)),j) |
112 | | -components(f::Fun{<:ContinuousSpace}) = components(Fun(f,canonicalspace(space(f)))) |
113 | | - |
114 | | - |
115 | | -function points(f::Fun{<:ContinuousSpace}) |
116 | | - n=ncoefficients(f) |
117 | | - d=domain(f) |
118 | | - K=ncomponents(d) |
119 | | - |
120 | | - m=isperiodic(d) ? max(K,n+2K-1) : n+K |
121 | | - points(f.space,m) |
122 | | -end |
123 | | - |
124 | | -## Conversion |
125 | | - |
126 | | -coefficients(cfsin::AbstractVector,A::ContinuousSpace,B::PiecewiseSpace) = |
127 | | - defaultcoefficients(cfsin,A,B) |
128 | | - |
129 | | -coefficients(cfsin::AbstractVector,A::PiecewiseSpace,B::ContinuousSpace) = |
130 | | - default_Fun(Fun(A,cfsin),B).coefficients |
131 | | - |
132 | | -coefficients(cfsin::AbstractVector,A::ContinuousSpace,B::ContinuousSpace) = |
133 | | - default_Fun(Fun(A,cfsin),B).coefficients |
134 | | - |
135 | 12 |
|
136 | 13 | # We implemnt conversion between continuous space and PiecewiseSpace with Chebyshev dirichlet |
137 | 14 |
|
@@ -407,34 +284,5 @@ function BlockBandedMatrix(S::SubOperator{T,<:ConcreteDirichlet{<:TensorChebyshe |
407 | 284 | end |
408 | 285 |
|
409 | 286 |
|
410 | | -union_rule(A::PiecewiseSpace, B::ContinuousSpace) = union(A, strictconvert(PiecewiseSpace, B)) |
411 | | -union_rule(A::ConstantSpace, B::ContinuousSpace) = B |
412 | 287 | union_rule(A::ContinuousSpace, B::PolynomialSpace{<:IntervalOrSegment}) = |
413 | 288 | Space(domain(A) ∪ domain(B)) |
414 | | - |
415 | | -function approx_union(a::AbstractVector{T}, b::AbstractVector{V}) where {T,V} |
416 | | - ret = sort!(union(a,b)) |
417 | | - for k in length(ret)-1:-1:1 |
418 | | - isapprox(ret[k] , ret[k+1]; atol=10eps()) && deleteat!(ret, k+1) |
419 | | - end |
420 | | - ret |
421 | | -end |
422 | | - |
423 | | - |
424 | | -function union_rule(A::ContinuousSpace{<:Real}, B::ContinuousSpace{<:Real}) |
425 | | - p_A,p_B = domain(A).points, domain(B).points |
426 | | - a,b = minimum(p_A), maximum(p_A) |
427 | | - c,d = minimum(p_B), maximum(p_B) |
428 | | - @assert !isempty((a..b) ∩ (c..d)) |
429 | | - ContinuousSpace(PiecewiseSegment(approx_union(p_A, p_B))) |
430 | | -end |
431 | | - |
432 | | -function integrate(f::Fun{<:ContinuousSpace}) |
433 | | - cs = [cumsum(x) for x in components(f)] |
434 | | - for k=1:length(cs)-1 |
435 | | - cs[k+1] += last(cs[k]) |
436 | | - end |
437 | | - Fun(Fun(cs, PiecewiseSpace), space(f)) |
438 | | -end |
439 | | - |
440 | | -cumsum(f::Fun{<:ContinuousSpace}) = integrate(f) |
0 commit comments