Mercurial > repos > public > sbplib_julia

changeset 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 79d350003339
children 1eefaefdd0c7
files DiffOps/src/laplace.jl SbpOperators/src/constantstenciloperator.jl
diffstat 2 files changed, 28 insertions(+), 54 deletions(-) [+]
line wrap: on
line diff
--- a/DiffOps/src/laplace.jl	Tue Jan 07 23:56:58 2020 +0100
+++ b/DiffOps/src/laplace.jl	Wed Jan 08 00:00:29 2020 +0100
@@ -42,7 +42,7 @@
 
 Implements the quadrature operator `H` of Dim dimension as a TensorMapping
 """
-struct Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
+struct Quadrature{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
     op::D2{T,N,M,K}
     grid::EquidistantGrid{Dim,T}
 end
@@ -51,8 +51,7 @@
 LazyTensors.range_size(H::Quadrature{2}, domain_size::NTuple{2,Integer}) where T = size(H.grid)
 LazyTensors.domain_size(H::Quadrature{2}, range_size::NTuple{2,Integer}) where T = size(H.grid)
 
-# TODO: Dispatch on Tuple{Index{R1},Index{R2}}?
-@inline function LazyTensors.apply(H::Quadrature{2}, v::AbstractArray{T,2} where T, I::Tuple{Index{R1}, Index{R2}}) where {R1, R2}
+@inline function LazyTensors.apply(H::Quadrature{2}, v::AbstractArray{T,2}, I::NTuple{2, Index}) where T
     N = size(H.grid)
     # Quadrature in x direction
     @inbounds q = apply_quadrature(H.op, spacing(H.grid)[1], v[I] , I[1], N[1])
@@ -61,24 +60,15 @@
     return q
 end
 
-function LazyTensors.apply(H::Quadrature{2}, v::AbstractArray{T,2} where T, i::NTuple{2,Integer})
-    I = Index{Unknown}.(i)
-    LazyTensors.apply(H, v, I)
-end
+LazyTensors.apply_transpose(H::Quadrature{2}, v::AbstractArray{T,2}, I::NTuple{2, Index}) where T = LazyTensors.apply(H,v,I)
 
-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)
-
-function LazyTensors.apply_transpose(H::Quadrature{2}, v::AbstractArray{T,2} where T, i::NTuple{2,Integer})
-    I = Index{Unknown}.(i)
-    LazyTensors.apply_transpose(H, v, I)
-end
 
 """
     InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
 
 Implements the inverse quadrature operator `inv(H)` of Dim dimension as a TensorMapping
 """
-struct InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
+struct InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
     op::D2{T,N,M,K}
     grid::EquidistantGrid{Dim,T}
 end
@@ -87,9 +77,7 @@
 LazyTensors.range_size(H_inv::InverseQuadrature{2}, domain_size::NTuple{2,Integer}) where T = size(H_inv.grid)
 LazyTensors.domain_size(H_inv::InverseQuadrature{2}, range_size::NTuple{2,Integer}) where T = size(H_inv.grid)
 
-# TODO: Dispatch on Tuple{Index{R1},Index{R2}}?
-@inline function LazyTensors.apply(H_inv::InverseQuadrature{2}, v::AbstractArray{T,2} where T, I::NTuple{2,Integer})
-    I = CartesianIndex(I);
+@inline function LazyTensors.apply(H_inv::InverseQuadrature{2}, v::AbstractArray{T,2}, I::NTuple{2, Index}) where T
     N = size(H_inv.grid)
     # Inverse quadrature in x direction
     @inbounds q_inv = apply_inverse_quadrature(H_inv.op, inverse_spacing(H_inv.grid)[1], v[I] , I[1], N[1])
@@ -98,7 +86,7 @@
     return q_inv
 end
 
-LazyTensors.apply_transpose(H_inv::InverseQuadrature{2}, v::AbstractArray{T,2} where T, I::NTuple{2,Integer}) = LazyTensors.apply(H_inv,v,I)
+LazyTensors.apply_transpose(H_inv::InverseQuadrature{2}, v::AbstractArray{T,2}, I::NTuple{2, Index}) where T = LazyTensors.apply(H_inv,v,I)
 
 """
     BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
@@ -118,15 +106,15 @@
 LazyTensors.domain_size(e::BoundaryValue{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
 
 # TODO: Make this independent of dimension
-function LazyTensors.apply(e::BoundaryValue, v::AbstractArray, I::NTuple{2,Int})
+function LazyTensors.apply(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{2, Index}) where T
     i = I[dim(e.bId)]
     j = I[3-dim(e.bId)]
     N_i = size(e.grid)[dim(e.bId)]
-    return apply_boundary_value(e.op, v[j], N_i, i, region(e.bId))
+    return apply_boundary_value(e.op, v[j], i, N_i, region(e.bId))
 end
 
-function LazyTensors.apply_transpose(e::BoundaryValue, v::AbstractArray, I::NTuple{1,Int})
-    u = selectdim(v,3-dim(e.bId),I[1])
+function LazyTensors.apply_transpose(e::BoundaryValue{T}, v::AbstractArray{T}, I::NTuple{1, Index}) where T
+    u = selectdim(v,3-dim(e.bId),Int(I[1]))
     return apply_boundary_value_transpose(e.op, u, region(e.bId))
 end
 
@@ -149,16 +137,16 @@
 
 # TODO: Not type stable D:<
 # TODO: Make this independent of dimension
-function LazyTensors.apply(d::NormalDerivative, v::AbstractArray, I::NTuple{2,Int})
+function LazyTensors.apply(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{2, Index}) where T
     i = I[dim(d.bId)]
     j = I[3-dim(d.bId)]
     N_i = size(d.grid)[dim(d.bId)]
     h_inv = inverse_spacing(d.grid)[dim(d.bId)]
-    return apply_normal_derivative(d.op, h_inv, v[j], N_i, i, region(d.bId))
+    return apply_normal_derivative(d.op, h_inv, v[j], i, N_i, region(d.bId))
 end
 
-function LazyTensors.apply_transpose(d::NormalDerivative, v::AbstractArray, I::NTuple{1,Int})
-    u = selectdim(v,3-dim(d.bId),I[1])
+function LazyTensors.apply_transpose(d::NormalDerivative{T}, v::AbstractArray{T}, I::NTuple{1, Index}) where T
+    u = selectdim(v,3-dim(d.bId),Int(I[1]))
     return apply_normal_derivative_transpose(d.op, inverse_spacing(d.grid)[dim(d.bId)], u, region(d.bId))
 end
 
@@ -176,13 +164,13 @@
 
 # TODO: Make this independent of dimension
 # TODO: Dispatch directly on Index{R}?
-function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Int}) where T
+function LazyTensors.apply(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1, Index}) where T
     h = spacing(q.grid)[3-dim(q.bId)]
     N = size(v)
     return apply_quadrature(q.op, h, v[I[1]], I[1], N[1])
 end
 
-LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1,Int}) where T = LazyTensors.apply(q,v,I)
+LazyTensors.apply_transpose(q::BoundaryQuadrature{T}, v::AbstractArray{T,1}, I::NTuple{1, Index}) where T = LazyTensors.apply(q,v,I)
 
 
 
--- a/SbpOperators/src/constantstenciloperator.jl	Tue Jan 07 23:56:58 2020 +0100
+++ b/SbpOperators/src/constantstenciloperator.jl	Wed Jan 08 00:00:29 2020 +0100
@@ -20,11 +20,6 @@
     i = Index(Int(index), r)
     return apply_2nd_derivative(op, h_inv, v, i)
 end
-
-# Wrapper functions for using regular indecies without specifying regions
-@inline function apply_2nd_derivative(op::ConstantStencilOperator, h_inv::Real, v::AbstractVector, i::Int)
-    return apply_2nd_derivative(op, h_inv, v, Index{Unknown}(i))
-end
 export apply_2nd_derivative
 
 apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Index{Lower}, N::Integer) where T = v*h*op.quadratureClosure[Int(i)]
@@ -36,11 +31,6 @@
     i = Index(Int(index), r)
     return apply_quadrature(op, h, v, i, N)
 end
-
-# Wrapper functions for using regular indecies without specifying regions
-function apply_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer) where T
-    return apply_quadrature(op, h, v, Index{Unknown}(i), N)
-end
 export apply_quadrature
 
 # TODO: Evaluate if divisions affect performance
@@ -54,10 +44,6 @@
     return apply_inverse_quadrature(op, h_inv, v, i, N)
 end
 
-# Wrapper functions for using regular indecies without specifying regions
-function apply_inverse_quadrature(op::ConstantStencilOperator, h::Real, v::T, i::Integer, N::Integer) where T
-    return apply_inverse_quadrature(op, h, v, Index{Unknown}(i), N)
-end
 export apply_inverse_quadrature
 
 function apply_boundary_value_transpose(op::ConstantStencilOperator, v::AbstractVector, ::Type{Lower})
@@ -75,18 +61,18 @@
 end
 export apply_boundary_value_transpose
 
-function apply_boundary_value(op::ConstantStencilOperator, v::Number, N::Integer, i::Integer, ::Type{Lower})
-    @boundscheck if !(0<length(i) <= N)
+function apply_boundary_value(op::ConstantStencilOperator, v::Number, i::Index, N::Integer, ::Type{Lower})
+    @boundscheck if !(0<length(Int(i)) <= N)
         throw(BoundsError())
     end
-    op.eClosure[i-1]*v
+    op.eClosure[Int(i)-1]*v
 end
 
-function apply_boundary_value(op::ConstantStencilOperator, v::Number, N::Integer, i::Integer, ::Type{Upper})
-    @boundscheck if !(0<length(i) <= N)
+function apply_boundary_value(op::ConstantStencilOperator, v::Number, i::Index,  N::Integer, ::Type{Upper})
+    @boundscheck if !(0<length(Int(i)) <= N)
         throw(BoundsError())
     end
-    op.eClosure[N-i]*v
+    op.eClosure[N-Int(i)]*v
 end
 export apply_boundary_value
 
@@ -106,18 +92,18 @@
 
 export apply_normal_derivative_transpose
 
-function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Lower})
-    @boundscheck if !(0<length(i) <= N)
+function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, i::Index, N::Integer, ::Type{Lower})
+    @boundscheck if !(0<length(Int(i)) <= N)
         throw(BoundsError())
     end
-    h_inv*op.dClosure[i-1]*v
+    h_inv*op.dClosure[Int(i)-1]*v
 end
 
-function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, N::Integer, i::Integer, ::Type{Upper})
-    @boundscheck if !(0<length(i) <= N)
+function apply_normal_derivative(op::ConstantStencilOperator, h_inv::Real, v::Number, i::Index, N::Integer, ::Type{Upper})
+    @boundscheck if !(0<length(Int(i)) <= N)
         throw(BoundsError())
     end
-    -h_inv*op.dClosure[N-i]*v
+    -h_inv*op.dClosure[N-Int(i)]*v
 end
 
 export apply_normal_derivative