Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/laplace/laplace.jl @ 333:01b851161018 refactor/combine_to_one_package

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Start converting to one package by moving all the files to their correct location
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 25 Sep 2020 13:06:02 +0200
parents SbpOperators/src/laplace/laplace.jl@9cc5d1498b2d
children 0844069ab5ff
comparison
equal deleted inserted replaced
332:535f1bff4bcc 333:01b851161018
1 export Laplace
2 """
3 Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
4
5 Implements the Laplace operator `L` in Dim dimensions as a tensor operator
6 The multi-dimensional tensor operator consists of a tuple of 1D SecondDerivative
7 tensor operators.
8 """
9 #export quadrature, inverse_quadrature, boundary_quadrature, boundary_value, normal_derivative
10 struct Laplace{Dim,T,N,M,K} <: TensorOperator{T,Dim}
11 D2::NTuple{Dim,SecondDerivative{T,N,M,K}}
12 #TODO: Write a good constructor
13 end
14
15 LazyTensors.domain_size(L::Laplace{Dim}, range_size::NTuple{Dim,Integer}) where {Dim} = range_size
16
17 function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim}
18 error("not implemented")
19 end
20
21 # u = L*v
22 function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, I::Index) where T
23 @inbounds u = LazyTensors.apply(L.D2[1],v,I)
24 return u
25 end
26
27 function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::Index, J::Index) where T
28 # 2nd x-derivative
29 @inbounds vx = view(v, :, Int(J))
30 @inbounds uᵢ = LazyTensors.apply(L.D2[1], vx , I)
31
32 # 2nd y-derivative
33 @inbounds vy = view(v, Int(I), :)
34 @inbounds uᵢ += LazyTensors.apply(L.D2[2], vy , J)
35
36 return uᵢ
37 end
38
39 LazyTensors.apply_transpose(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::Vararg{Index,Dim}) where {T,Dim} = LazyTensors.apply(L, v, I...)
40
41 # quadrature(L::Laplace) = Quadrature(L.op, L.grid)
42 # inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid)
43 # boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId)
44 # normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId)
45 # boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId)
46 # export NormalDerivative
47 # """
48 # NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
49 #
50 # Implements the boundary operator `d` as a TensorMapping
51 # """
52 # struct NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
53 # op::D2{T,N,M,K}
54 # grid::EquidistantGrid{2}
55 # bId::CartesianBoundary
56 # end
57 #
58 # # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue?
59 # # Can we give special treatment to TensorMappings that go to a higher dim?
60 # function LazyTensors.range_size(e::NormalDerivative, domain_size::NTuple{1,Integer})
61 # if dim(e.bId) == 1
62 # return (UnknownDim, domain_size[1])
63 # elseif dim(e.bId) == 2
64 # return (domain_size[1], UnknownDim)
65 # end
66 # end
67 # LazyTensors.domain_size(e::NormalDerivative, range_size::NTuple{2,Integer}) = (range_size[3-dim(e.bId)],)
68 #
69 # # TODO: Not type stable D:<
70 # # TODO: Make this independent of dimension
71 # function LazyTensors.apply(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{2,Index}) where T
72 # i = I[dim(d.bId)]
73 # j = I[3-dim(d.bId)]
74 # N_i = size(d.grid)[dim(d.bId)]
75 # h_inv = inverse_spacing(d.grid)[dim(d.bId)]
76 # return apply_normal_derivative(d.op, h_inv, v[j], i, N_i, region(d.bId))
77 # end
78 #
79 # function LazyTensors.apply_transpose(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{1,Index}) where T
80 # u = selectdim(v,3-dim(d.bId),Int(I[1]))
81 # return apply_normal_derivative_transpose(d.op, inverse_spacing(d.grid)[dim(d.bId)], u, region(d.bId))
82 # end
83 #
84 # """
85 # BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1}
86 #
87 # Implements the boundary operator `q` as a TensorOperator
88 # """
89 # export BoundaryQuadrature
90 # struct BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1}
91 # op::D2{T,N,M,K}
92 # grid::EquidistantGrid{2}
93 # bId::CartesianBoundary
94 # end
95 #
96 #
97 # # TODO: Make this independent of dimension
98 # function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T
99 # h = spacing(q.grid)[3-dim(q.bId)]
100 # N = size(v)
101 # return apply_quadrature(q.op, h, v[I[1]], I[1], N[1])
102 # end
103 #
104 # LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Index}) where T = LazyTensors.apply(q,v,I)
105 #
106 #
107 #
108 #
109 # struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end
110 #
111 # function sat(L::Laplace{2,T}, bc::Neumann{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, I::CartesianIndex{2}) where {T,Bid}
112 # e = boundary_value(L, Bid())
113 # d = normal_derivative(L, Bid())
114 # Hᵧ = boundary_quadrature(L, Bid())
115 # H⁻¹ = inverse_quadrature(L)
116 # return (-H⁻¹*e*Hᵧ*(d'*v - g))[I]
117 # end
118 #
119 # struct Dirichlet{Bid<:BoundaryIdentifier} <: BoundaryCondition
120 # tau::Float64
121 # end
122 #
123 # function sat(L::Laplace{2,T}, bc::Dirichlet{Bid}, v::AbstractArray{T,2}, g::AbstractVector{T}, i::CartesianIndex{2}) where {T,Bid}
124 # e = boundary_value(L, Bid())
125 # d = normal_derivative(L, Bid())
126 # Hᵧ = boundary_quadrature(L, Bid())
127 # H⁻¹ = inverse_quadrature(L)
128 # return (-H⁻¹*(tau/h*e + d)*Hᵧ*(e'*v - g))[I]
129 # # Need to handle scalar multiplication and addition of TensorMapping
130 # end
131
132 # function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D
133 # return apply(s.L, v, i) +
134 # sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) +
135 # sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) +
136 # sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) +
137 # sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i)
138 # end