Mercurial > repos > public > sbplib_julia

comparison SbpOperators/src/laplace/laplace.jl @ 304:5645021683d3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 09 Sep 2020 20:41:31 +0200
parents 6fa2ba769ae3
children c1fcc35e19cb
comparison
equal deleted inserted replaced
303:d5475ad78b28 304:5645021683d3
1 """ 1 """
2 Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim} 2 Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
3 3
4 Implements the Laplace operator `L` in Dim dimensions as a tensor operator 4 Implements the Laplace operator `L` in Dim dimensions as a tensor operator
5 The multi-dimensional tensor operator simply consists of a tuple of the 1D 5 The multi-dimensional tensor operator consists of a tuple of 1D SecondDerivative
6 Laplace tensor operator as defined by ConstantLaplaceOperator. 6 tensor operators.
7 """ 7 """
8 struct Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim} 8 struct Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
9 D2::NTuple(Dim,SecondDerivative{T,N,M,K}) 9 D2::NTuple(Dim,SecondDerivative{T,N,M,K})
10 #TODO: Write a good constructor 10 #TODO: Write a good constructor
11 end 11 end
15 15
16 function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim} 16 function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim}
17 error("not implemented") 17 error("not implemented")
18 end 18 end
19 19
20 function LazyTensors.apply_transpose(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim} = LazyTensors.apply(L, v, I)
21
20 # u = L*v 22 # u = L*v
21 function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T 23 function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
22 return apply(L.D2[1],v,I) 24 @inbounds u = apply(L.D2[1],v,I)
25 return u
23 end 26 end
24 27
25 28
26 @inline function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T 29 @inline function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
27 # 2nd x-derivative 30 # 2nd x-derivative
28 @inbounds vx = view(v, :, Int(I[2])) 31 @inbounds vx = view(v, :, Int(I[2]))
29 @inbounds uᵢ = apply(L.D2[1], vx , (I[1],)) #Tuple conversion here is ugly. How to do it? Should we use indexing here? 32 @inbounds uᵢ = apply(L.D2[1], vx , I[1])
30 33
31 # 2nd y-derivative 34 # 2nd y-derivative
32 @inbounds vy = view(v, Int(I[1]), :) 35 @inbounds vy = view(v, Int(I[1]), :)
33 @inbounds uᵢ += apply(L.D2[2], vy , (I[2],)) #Tuple conversion here is ugly. How to do it? 36 @inbounds uᵢ += apply(L.D2[2], vy , I[2])
34 37
35 return uᵢ 38 return uᵢ
36 end 39 end
37 40
38 quadrature(L::Laplace) = Quadrature(L.op, L.grid) 41 quadrature(L::Laplace) = Quadrature(L.op, L.grid)
39 inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid) 42 inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid)
40 boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId) 43 boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId)
41 normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId) 44 normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId)
42 boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId) 45 boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId)
43 export quadrature 46 export quadrature
44
45 # At the moment the grid property is used all over. It could possibly be removed if we implement all the 1D operators as TensorMappings
46 """
47 Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
48
49 Implements the quadrature operator `H` of Dim dimension as a TensorMapping
50 """
51 struct Quadrature{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
52 op::D2{T,N,M,K}
53 grid::EquidistantGrid{Dim,T}
54 end
55 export Quadrature
56
57 LazyTensors.domain_size(H::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
58
59 @inline function LazyTensors.apply(H::Quadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
60 N = size(H.grid)
61 # Quadrature in x direction
62 @inbounds q = apply_quadrature(H.op, spacing(H.grid)[1], v[I] , I[1], N[1])
63 # Quadrature in y-direction
64 @inbounds q = apply_quadrature(H.op, spacing(H.grid)[2], q, I[2], N[2])
65 return q
66 end
67
68 LazyTensors.apply_transpose(H::Quadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T = LazyTensors.apply(H,v,I)
69
70
71 """
72 InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
73
74 Implements the inverse quadrature operator `inv(H)` of Dim dimension as a TensorMapping
75 """
76 struct InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
77 op::D2{T,N,M,K}
78 grid::EquidistantGrid{Dim,T}
79 end
80 export InverseQuadrature
81
82 LazyTensors.domain_size(H_inv::InverseQuadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
83
84 @inline function LazyTensors.apply(H_inv::InverseQuadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
85 N = size(H_inv.grid)
86 # Inverse quadrature in x direction
87 @inbounds q_inv = apply_inverse_quadrature(H_inv.op, inverse_spacing(H_inv.grid)[1], v[I] , I[1], N[1])
88 # Inverse quadrature in y-direction
89 @inbounds q_inv = apply_inverse_quadrature(H_inv.op, inverse_spacing(H_inv.grid)[2], q_inv, I[2], N[2])
90 return q_inv
91 end
92
93 LazyTensors.apply_transpose(H_inv::InverseQuadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T = LazyTensors.apply(H_inv,v,I)
94 47
95 """ 48 """
96 BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1} 49 BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
97 50
98 Implements the boundary operator `e` as a TensorMapping 51 Implements the boundary operator `e` as a TensorMapping