Mercurial > repos > public > sbplib_julia

comparison DiffOps/src/laplace.jl @ 280:fe9e8737ddfa boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Change to using region indices in apply of BoundaryValue, NormalDerivative and BoundaryQuadrature
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 08 Jan 2020 00:00:29 +0100
parents babc4288e6a6
children 1eefaefdd0c7
comparison
equal deleted inserted replaced
279:79d350003339 280:fe9e8737ddfa
40 """ 40 """
41 Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} 41 Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
42 42
43 Implements the quadrature operator `H` of Dim dimension as a TensorMapping 43 Implements the quadrature operator `H` of Dim dimension as a TensorMapping
44 """ 44 """
45 struct Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} 45 struct Quadrature{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
46 op::D2{T,N,M,K} 46 op::D2{T,N,M,K}
47 grid::EquidistantGrid{Dim,T} 47 grid::EquidistantGrid{Dim,T}
48 end 48 end
49 export Quadrature 49 export Quadrature
50 50
51 LazyTensors.range_size(H::Quadrature{2}, domain_size::NTuple{2,Integer}) where T = size(H.grid) 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) 52 LazyTensors.domain_size(H::Quadrature{2}, range_size::NTuple{2,Integer}) where T = size(H.grid)
53 53
54 # TODO: Dispatch on Tuple{Index{R1},Index{R2}}? 54 @inline function LazyTensors.apply(H::Quadrature{2}, v::AbstractArray{T,2}, I::NTuple{2, Index}) where T
55 @inline function LazyTensors.apply(H::Quadrature{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
56 N = size(H.grid) 55 N = size(H.grid)
57 # Quadrature in x direction 56 # Quadrature in x direction
58 @inbounds q = apply_quadrature(H.op, spacing(H.grid)[1], v[I] , I[1], N[1]) 57 @inbounds q = apply_quadrature(H.op, spacing(H.grid)[1], v[I] , I[1], N[1])
59 # Quadrature in y-direction 58 # Quadrature in y-direction
60 @inbounds q = apply_quadrature(H.op, spacing(H.grid)[2], q, I[2], N[2]) 59 @inbounds q = apply_quadrature(H.op, spacing(H.grid)[2], q, I[2], N[2])
61 return q 60 return q
62 end 61 end
63 62
64 function LazyTensors.apply(H::Quadrature{2}, v::AbstractArray{T,2} where T, i::NTuple{2,Integer}) 63 LazyTensors.apply_transpose(H::Quadrature{2}, v::AbstractArray{T,2}, I::NTuple{2, Index}) where T = LazyTensors.apply(H,v,I)
65 I = Index{Unknown}.(i) 64
66 LazyTensors.apply(H, v, I)
67 end
68
69 LazyTensors.apply_transpose(H::Quadrature{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}} where {R1, R2}) = LazyTensors.apply(H,v,I)
70
71 function LazyTensors.apply_transpose(H::Quadrature{2}, v::AbstractArray{T,2} where T, i::NTuple{2,Integer})
72 I = Index{Unknown}.(i)
73 LazyTensors.apply_transpose(H, v, I)
74 end
75 65
76 """ 66 """
77 InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} 67 InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
78 68
79 Implements the inverse quadrature operator `inv(H)` of Dim dimension as a TensorMapping 69 Implements the inverse quadrature operator `inv(H)` of Dim dimension as a TensorMapping
80 """ 70 """
81 struct InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} 71 struct InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
82 op::D2{T,N,M,K} 72 op::D2{T,N,M,K}
83 grid::EquidistantGrid{Dim,T} 73 grid::EquidistantGrid{Dim,T}
84 end 74 end
85 export InverseQuadrature 75 export InverseQuadrature
86 76
87 LazyTensors.range_size(H_inv::InverseQuadrature{2}, domain_size::NTuple{2,Integer}) where T = size(H_inv.grid) 77 LazyTensors.range_size(H_inv::InverseQuadrature{2}, domain_size::NTuple{2,Integer}) where T = size(H_inv.grid)
88 LazyTensors.domain_size(H_inv::InverseQuadrature{2}, range_size::NTuple{2,Integer}) where T = size(H_inv.grid) 78 LazyTensors.domain_size(H_inv::InverseQuadrature{2}, range_size::NTuple{2,Integer}) where T = size(H_inv.grid)
89 79
90 # TODO: Dispatch on Tuple{Index{R1},Index{R2}}? 80 @inline function LazyTensors.apply(H_inv::InverseQuadrature{2}, v::AbstractArray{T,2}, I::NTuple{2, Index}) where T
91 @inline function LazyTensors.apply(H_inv::InverseQuadrature{2}, v::AbstractArray{T,2} where T, I::NTuple{2,Integer})
92 I = CartesianIndex(I);
93 N = size(H_inv.grid) 81 N = size(H_inv.grid)
94 # Inverse quadrature in x direction 82 # Inverse quadrature in x direction
95 @inbounds q_inv = apply_inverse_quadrature(H_inv.op, inverse_spacing(H_inv.grid)[1], v[I] , I[1], N[1]) 83 @inbounds q_inv = apply_inverse_quadrature(H_inv.op, inverse_spacing(H_inv.grid)[1], v[I] , I[1], N[1])
96 # Inverse quadrature in y-direction 84 # Inverse quadrature in y-direction
97 @inbounds q_inv = apply_inverse_quadrature(H_inv.op, inverse_spacing(H_inv.grid)[2], q_inv, I[2], N[2]) 85 @inbounds q_inv = apply_inverse_quadrature(H_inv.op, inverse_spacing(H_inv.grid)[2], q_inv, I[2], N[2])
98 return q_inv 86 return q_inv
99 end 87 end
100 88
101 LazyTensors.apply_transpose(H_inv::InverseQuadrature{2}, v::AbstractArray{T,2} where T, I::NTuple{2,Integer}) = LazyTensors.apply(H_inv,v,I) 89 LazyTensors.apply_transpose(H_inv::InverseQuadrature{2}, v::AbstractArray{T,2}, I::NTuple{2, Index}) where T = LazyTensors.apply(H_inv,v,I)
102 90
103 """ 91 """
104 BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1} 92 BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
105 93
106 Implements the boundary operator `e` as a TensorMapping 94 Implements the boundary operator `e` as a TensorMapping
116 # Can we give special treatment to TensorMappings that go to a higher dim? 104 # Can we give special treatment to TensorMappings that go to a higher dim?
117 LazyTensors.range_size(e::BoundaryValue{T}, domain_size::NTuple{1,Integer}) where T = size(e.grid) 105 LazyTensors.range_size(e::BoundaryValue{T}, domain_size::NTuple{1,Integer}) where T = size(e.grid)
118 LazyTensors.domain_size(e::BoundaryValue{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],) 106 LazyTensors.domain_size(e::BoundaryValue{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
119 107
120 # TODO: Make this independent of dimension 108 # TODO: Make this independent of dimension
121 function LazyTensors.apply(e::BoundaryValue, v::AbstractArray, I::NTuple{2,Int}) 109 function LazyTensors.apply(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{2, Index}) where T
122 i = I[dim(e.bId)] 110 i = I[dim(e.bId)]
123 j = I[3-dim(e.bId)] 111 j = I[3-dim(e.bId)]
124 N_i = size(e.grid)[dim(e.bId)] 112 N_i = size(e.grid)[dim(e.bId)]
125 return apply_boundary_value(e.op, v[j], N_i, i, region(e.bId)) 113 return apply_boundary_value(e.op, v[j], i, N_i, region(e.bId))
126 end 114 end
127 115
128 function LazyTensors.apply_transpose(e::BoundaryValue, v::AbstractArray, I::NTuple{1,Int}) 116 function LazyTensors.apply_transpose(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{1, Index}) where T
129 u = selectdim(v,3-dim(e.bId),I[1]) 117 u = selectdim(v,3-dim(e.bId),Int(I[1]))
130 return apply_boundary_value_transpose(e.op, u, region(e.bId)) 118 return apply_boundary_value_transpose(e.op, u, region(e.bId))
131 end 119 end
132 120
133 """ 121 """
134 NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1} 122 NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
147 LazyTensors.range_size(e::NormalDerivative{T}, domain_size::NTuple{1,Integer}) where T = size(e.grid) 135 LazyTensors.range_size(e::NormalDerivative{T}, domain_size::NTuple{1,Integer}) where T = size(e.grid)
148 LazyTensors.domain_size(e::NormalDerivative{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],) 136 LazyTensors.domain_size(e::NormalDerivative{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
149 137
150 # TODO: Not type stable D:< 138 # TODO: Not type stable D:<
151 # TODO: Make this independent of dimension 139 # TODO: Make this independent of dimension
152 function LazyTensors.apply(d::NormalDerivative, v::AbstractArray, I::NTuple{2,Int}) 140 function LazyTensors.apply(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{2, Index}) where T
153 i = I[dim(d.bId)] 141 i = I[dim(d.bId)]
154 j = I[3-dim(d.bId)] 142 j = I[3-dim(d.bId)]
155 N_i = size(d.grid)[dim(d.bId)] 143 N_i = size(d.grid)[dim(d.bId)]
156 h_inv = inverse_spacing(d.grid)[dim(d.bId)] 144 h_inv = inverse_spacing(d.grid)[dim(d.bId)]
157 return apply_normal_derivative(d.op, h_inv, v[j], N_i, i, region(d.bId)) 145 return apply_normal_derivative(d.op, h_inv, v[j], i, N_i, region(d.bId))
158 end 146 end
159 147
160 function LazyTensors.apply_transpose(d::NormalDerivative, v::AbstractArray, I::NTuple{1,Int}) 148 function LazyTensors.apply_transpose(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{1, Index}) where T
161 u = selectdim(v,3-dim(d.bId),I[1]) 149 u = selectdim(v,3-dim(d.bId),Int(I[1]))
162 return apply_normal_derivative_transpose(d.op, inverse_spacing(d.grid)[dim(d.bId)], u, region(d.bId)) 150 return apply_normal_derivative_transpose(d.op, inverse_spacing(d.grid)[dim(d.bId)], u, region(d.bId))
163 end 151 end
164 152
165 """ 153 """
166 BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1} 154 BoundaryQuadrature{T,N,M,K} <: TensorOperator{T,1}
174 end 162 end
175 export BoundaryQuadrature 163 export BoundaryQuadrature
176 164
177 # TODO: Make this independent of dimension 165 # TODO: Make this independent of dimension
178 # TODO: Dispatch directly on Index{R}? 166 # TODO: Dispatch directly on Index{R}?
179 function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Int}) where T 167 function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1, Index}) where T
180 h = spacing(q.grid)[3-dim(q.bId)] 168 h = spacing(q.grid)[3-dim(q.bId)]
181 N = size(v) 169 N = size(v)
182 return apply_quadrature(q.op, h, v[I[1]], I[1], N[1]) 170 return apply_quadrature(q.op, h, v[I[1]], I[1], N[1])
183 end 171 end
184 172
185 LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Int}) where T = LazyTensors.apply(q,v,I) 173 LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1, Index}) where T = LazyTensors.apply(q,v,I)
186 174
187 175
188 176
189 177
190 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end 178 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end