Mercurial > repos > public > sbplib_julia

comparison DiffOps/src/laplace.jl @ 262:f1e90a92ad74 boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add Quadrature and InverseQuadrature for Laplace as TensorMappings. Implement and test the 2D case. Fix implementation of apply_transpose for BoundaryQuadrature and add tests.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 26 Nov 2019 08:28:26 -0800
parents 5571d2c5bf0f
children 9ad447176ba1
comparison
equal deleted inserted replaced
261:01017d2b46b0 262:f1e90a92ad74
29 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2}) 29 function apply(L::Laplace{2}, v::AbstractArray{T,2} where T, i::CartesianIndex{2})
30 I = Index{Unknown}.(Tuple(i)) 30 I = Index{Unknown}.(Tuple(i))
31 apply(L, v, I) 31 apply(L, v, I)
32 end 32 end
33 33
34 quadrature(L::Laplace) = Quadrature(L.op, L.grid)
35 inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid)
34 boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId) 36 boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId)
35 normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId) 37 normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId)
36 boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId) 38 boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId)
37 39
40 """
41 Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
42
43 Implements the quadrature operator `H` of Dim dimension as a TensorMapping
44 """
45 struct Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
46 op::D2{T,N,M,K}
47 grid::EquidistantGrid{Dim,T}
48 end
49 export Quadrature
50
51 LazyTensors.range_size(H::Quadrature{2}, domain_size::NTuple{2,Integer}) where T = size(H.grid)
52 LazyTensors.domain_size(H::Quadrature{2}, range_size::NTuple{2,Integer}) where T = size(H.grid)
53
54 # TODO: Dispatch on Tuple{Index{R1},Index{R2}}?
55 @inline function LazyTensors.apply(H::Quadrature{2}, v::AbstractArray{T,2} where T, I::NTuple{2,Integer})
56 I = CartesianIndex(I);
57 N = size(H.grid)
58 # Quadrature in x direction
59 @inbounds q = apply_quadrature(H.op, H.grid.spacing[1], v[I] , I[1], N[1])
60 # Quadrature in y-direction
61 @inbounds q = apply_quadrature(H.op, H.grid.spacing[2], q, I[2], N[2])
62 return q
63 end
64
65 LazyTensors.apply_transpose(H::Quadrature{2}, v::AbstractArray{T,2} where T, I::NTuple{2,Integer}) = LazyTensors.apply(H,v,I)
66
67 """
68 InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
69
70 Implements the inverse quadrature operator `inv(H)` of Dim dimension as a TensorMapping
71 """
72 struct InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
73 op::D2{T,N,M,K}
74 grid::EquidistantGrid{Dim,T}
75 end
76 export InverseQuadrature
77
78 LazyTensors.range_size(H_inv::InverseQuadrature{2}, domain_size::NTuple{2,Integer}) where T = size(H_inv.grid)
79 LazyTensors.domain_size(H_inv::InverseQuadrature{2}, range_size::NTuple{2,Integer}) where T = size(H_inv.grid)
80
81 # TODO: Dispatch on Tuple{Index{R1},Index{R2}}?
82 @inline function LazyTensors.apply(H_inv::InverseQuadrature{2}, v::AbstractArray{T,2} where T, I::NTuple{2,Integer})
83 I = CartesianIndex(I);
84 N = size(H_inv.grid)
85 # Inverse quadrature in x direction
86 @inbounds q_inv = apply_inverse_quadrature(H_inv.op, H_inv.grid.inverse_spacing[1], v[I] , I[1], N[1])
87 # Inverse quadrature in y-direction
88 @inbounds q_inv = apply_inverse_quadrature(H_inv.op, H_inv.grid.inverse_spacing[2], q_inv, I[2], N[2])
89 return q_inv
90 end
91
92 LazyTensors.apply_transpose(H_inv::InverseQuadrature{2}, v::AbstractArray{T,2} where T, I::NTuple{2,Integer}) = LazyTensors.apply(H_inv,v,I)
38 93
39 """ 94 """
40 BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1} 95 BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
41 96
42 Implements the boundary operator `e` as a TensorMapping 97 Implements the boundary operator `e` as a TensorMapping
63 118
64 function LazyTensors.apply_transpose(e::BoundaryValue, v::AbstractArray, I::NTuple{1,Int}) 119 function LazyTensors.apply_transpose(e::BoundaryValue, v::AbstractArray, I::NTuple{1,Int})
65 u = selectdim(v,3-dim(e.bId),I[1]) 120 u = selectdim(v,3-dim(e.bId),I[1])
66 return apply_e_T(e.op, u, region(e.bId)) 121 return apply_e_T(e.op, u, region(e.bId))
67 end 122 end
68
69
70 123
71 """ 124 """
72 NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1} 125 NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
73 126
74 Implements the boundary operator `d` as a TensorMapping 127 Implements the boundary operator `d` as a TensorMapping
111 bId::CartesianBoundary 164 bId::CartesianBoundary
112 end 165 end
113 export BoundaryQuadrature 166 export BoundaryQuadrature
114 167
115 # TODO: Make this independent of dimension 168 # TODO: Make this independent of dimension
169 # TODO: Dispatch directly on Index{R}?
116 function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Int}) where T 170 function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Int}) where T
117 h = spacing(q.grid)[3-dim(q.bId)] 171 h = q.grid.spacing[3-dim(q.bId)]
118 N = size(v) 172 N = size(v)
119 return apply_quadrature(q.op, h, v[I[1]], I[1], N[1]) 173 return apply_quadrature(q.op, h, v[I[1]], I[1], N[1])
120 end 174 end
121 175
122 LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Int}) where T = apply(q,v,I) 176 LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Int}) where T = LazyTensors.apply(q,v,I)
123 177
124 178
125 179
126 180
127 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end 181 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end