Mercurial > repos > public > sbplib_julia

comparison sbpD2.jl @ 43:ef060ab3b035

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
remove stride and remove some bugs
author Ylva Rydin <ylva.rydin@telia.com>
date Thu, 10 Jan 2019 16:30:54 +0100
parents 8b04efde1a46
children 8c4cfa680e5a
comparison
equal deleted inserted replaced
40:8b04efde1a46 43:ef060ab3b035
1 abstract type ConstantStencilOperator end 1 abstract type ConstantStencilOperator end
2 2
3 function apply!(op::ConstantStencilOperator, u::AbstractVector, v::AbstractVector, h::Real, start::Int, N::Int; stride::Int=1) 3 function apply!(op::ConstantStencilOperator, u::AbstractVector, v::AbstractVector, h::Real)
4 N = length(v)
4 cSize = closureSize(op) 5 cSize = closureSize(op)
5 6
6 for i ∈ range(start; length=cSize, step=stride) 7 for i ∈ range(1; length=cSize)
7 u[i] = apply!(op.closureStencils[i], v, i; stride=stride)/h^2 8 u[i] = apply(op.closureStencils[i], v, i)/h^2
8 end 9 end
9 10
10 innerStart = start + cSize*stride 11 innerStart = 1 + cSize
11 innerEnd = N - cSize*stride-1 12 innerEnd = N - cSize
12 for i ∈ range(innerStart, stop=innerEnd, step=stride) 13 for i ∈ range(innerStart, stop=innerEnd)
13 u[i] = apply(op.innerStencil, v, i; stride=stride)/h^2 14 u[i] = apply(op.innerStencil, v, i)/h^2
14 end 15 end
15 16
16 for i ∈ range(innerEnd+1, length=cSize, step=cSize) 17 for i ∈ range(innerEnd+1, length=cSize)
17 u[i] = op.parity*apply(flip(op.closureStencils[M-i+1]), v, i; stride=stride)/h^2 18 u[i] = op.parity*apply(flip(op.closureStencils[N-i+1]), v, i)/h^2
18 end 19 end
19 end 20 end
20 21
21 odd = -1 22 odd = -1
22 even = 1 23 even = 1
29 dClosure::Vector{T} 30 dClosure::Vector{T}
30 parity::Int 31 parity::Int
31 end 32 end
32 33
33 function closureSize(D::D2)::Int 34 function closureSize(D::D2)::Int
34 return length(quadratureClosure) 35 return length(D.quadratureClosure)
35 end 36 end
36 37
37 function readOperator(D2fn, Hfn) 38 function readOperator(D2fn, Hfn)
38 d = readSectionedFile(D2fn) 39 d = readSectionedFile(D2fn)
39 h = readSectionedFile(Hfn) 40 h = readSectionedFile(Hfn)