Mercurial > repos > public > sbplib_julia

comparison diffOp.jl @ 168:45840a8127d6 boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix typos and add NormalDerivative
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 21 May 2019 14:20:20 +0200
parents 89b63bdf1ea8
children 24ee4def7ffb
comparison
equal deleted inserted replaced
167:074cb898f24e 168:45840a8127d6
11 11
12 function matrixRepresentation(D::DiffOp) 12 function matrixRepresentation(D::DiffOp)
13 error("not implemented") 13 error("not implemented")
14 end 14 end
15 15
16 function boundaryCondition(D::DiffOp,b::Grid.BoundaryId,type)::(Closure, Penalty) 16 function boundaryCondition(D::DiffOp,b::BoundaryIdentifier,type)::(Closure, Penalty)
17 error("not implemented") 17 error("not implemented")
18 end 18 end
19 19
20 function interface(Du::DiffOp, Dv::DiffOp, b::Grid.BoundaryId; type) 20 function interface(Du::DiffOp, Dv::DiffOp, b::BoundaryIdentifier; type)
21 error("not implemented") 21 error("not implemented")
22 end 22 end
23 23
24 abstract type Closure end 24 abstract type Closure end
25 25
96 u = zeros(eltype(v), size(v)) 96 u = zeros(eltype(v), size(v))
97 apply!(D,v,u) 97 apply!(D,v,u)
98 return u 98 return u
99 end 99 end
100 100
101 struct NormalDerivative{N,M,K}
102 op::D2{Float64,N,M,K}
103 grid::EquidistantGrid
104 bId::CartesianBoundary
105 end
106
101 struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim} 107 struct Laplace{Dim,T<:Real,N,M,K} <: DiffOpCartesian{Dim}
102 grid::Grid.EquidistantGrid{Dim,T} 108 grid::EquidistantGrid{Dim,T}
103 a::T 109 a::T
104 op::D2{Float64,N,M,K} 110 op::D2{Float64,N,M,K}
105 e::BoundaryValue 111 # e::BoundaryValue
106 d::NormalDerivative 112 d::NormalDerivative
107 end 113 end
108 114
109 function apply(L::Laplace{Dim}, v::AbstractArray{T,Dim} where T, I::CartesianIndex{Dim}) where Dim 115 function apply(L::Laplace{Dim}, v::AbstractArray{T,Dim} where T, I::CartesianIndex{Dim}) where Dim
110 error("not implemented") 116 error("not implemented")
134 140
135 struct BoundaryOperator 141 struct BoundaryOperator
136 142
137 end 143 end
138 144
139 struct BoundaryValue 145 struct BoundaryValue{N,M,K}
140 op::D2{Float64,N,M,K} 146 op::D2{Float64,N,M,K}
141 end 147 end
142 148
143 function apply(e::BoundaryValue) 149 function apply(e::BoundaryValue)
144 150
146 152
147 function apply_adjoint(e::BoundaryValue) 153 function apply_adjoint(e::BoundaryValue)
148 154
149 end 155 end
150 156
151 struct NormalDerivative 157
152 op::D2{Float64,N,M,K} 158
153 end 159 function apply_transpose(d::NormalDerivative, v::AbstractArray, I::Integer)
154 160 u = selectdim(v,dim(d.bId),I)
155 function apply(e::NormalDerivative) 161 return apply_d(d.op, h, u, region(d.bId))
156 162 end
157 end 163
158 164 # Not correct abstraction level
159 function apply_adjoint(e::NormalDerivative) 165 # TODO: Not type stable D:<
160 166 function apply(d::NormalDerivative, v::AbstractArray, I::Tuple{Integer,Integer})
167 i = I[dim(d.bId)]
168 j = I[3-dim(d.bId)]
169 N_i = d.grid.size[dim(d.bId)]
170
171 r = getregion(i, closureSize(d.op), N_i)
172
173 if r != region(d.bId)
174 return 0
175 end
176
177 if r == Lower
178 # Note, closures are indexed by offset. Fix this D:<
179 return d.op.dClosure[i-1]*v[j]
180 elseif r == Upper
181 return d.op.dClosure[N_i-j]*v[j]
182 end
161 end 183 end
162 184
163 """ 185 """
164 A BoundaryCondition should implement the method 186 A BoundaryCondition should implement the method
165 sat(::DiffOp, v::AbstractArray, data::AbstractArray, ...) 187 sat(::DiffOp, v::AbstractArray, data::AbstractArray, ...)
171 end 193 end
172 194
173 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition 195 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition
174 end 196 end
175 197
176 function sat(L::Laplace{2}, bc::Neumann{CartesianBoundary{1,R}}, v::AbstractArray{T,2} where T, g::AbstractVector{T}, i::CartesianIndex{2}) where R 198 function sat(L::Laplace{2}, bc::Neumann{CartesianBoundary{1,R}}, v::AbstractArray{T,2} where T, g::AbstractVector, i::CartesianIndex{2}) where R
177 199
178 # Hi * e * H_gamma * (d'*v - g) 200 # Hi * e * H_gamma * (d'*v - g)
179 # e, d, H_gamma applied based on bc.boundaryId 201 # e, d, H_gamma applied based on bc.boundaryId
180 end 202 end
181 203
182 function sat(L::Laplace{2}, bc::Dirichlet{CartesianBoundary{1,R}}, v::AbstractArray{T,2} where T, g::AbstractVector{T}, i::CartesianIndex{2}) where R 204 function sat(L::Laplace{2}, bc::Dirichlet{CartesianBoundary{1,R}}, v::AbstractArray{T,2} where T, g::AbstractVector, i::CartesianIndex{2}) where R
183 # Hi * (tau/h*e + sig*d) * H_gamma * (e'*v - g) 205 # Hi * (tau/h*e + sig*d) * H_gamma * (e'*v - g)
184 # e, d, H_gamma applied based on bc.boundaryId 206 # e, d, H_gamma applied based on bc.boundaryId
185 end 207 end
186 208
187 function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D 209 # function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D
188 return apply(s.L, v, i) + 210 # return apply(s.L, v, i) +
189 sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) + 211 # sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) +
190 sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) + 212 # sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) +
191 sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) + 213 # sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) +
192 sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i) + 214 # sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i)
193 end 215 # end