Mercurial > repos > public > sbplib_julia

comparison DiffOps/src/laplace.jl @ 286:7247e85dc1e8 tensor_mappings

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Start separating ConstantStencilOp into multiple 1D tensor mappings, e.g. ConstantLaplaceOp. Sketch an implementation of the multi-D laplace tensor operator as a tuple of 1D laplace tensor operators.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 22 Jun 2020 19:47:20 +0200
parents e21dcda55163
children
comparison
equal deleted inserted replaced
285:e21dcda55163 286:7247e85dc1e8
1 #TODO: move to sbpoperators.jl
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 simply consists of a tuple of the 1D
7 Laplace tensor operator as defined by ConstantLaplaceOperator.
8 """
1 struct Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim} 9 struct Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
2 grid::EquidistantGrid{Dim,T} # TODO: Should this be here? Should probably be possible to applay a Laplace object to any size grid function 10 tensorOps::NTuple(Dim,ConstantLaplaceOperator{T,N,M,K})
3 a::T # TODO: Better name? 11 #TODO: Write a good constructor
4 op::D2{T,N,M,K}
5 end 12 end
6 export Laplace 13 export Laplace
7 14
8 # At the moment the grid property is used all over. It could possibly be removed if we implement all the 1D operators as TensorMappings 15 LazyTensors.domain_size(H::Laplace{Dim}, range_size::NTuple{Dim,Integer}) = range_size
9
10 LazyTensors.domain_size(H::Laplace{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
11 16
12 function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim} 17 function LazyTensors.apply(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim}
13 error("not implemented") 18 error("not implemented")
14 end 19 end
15 20
16 # u = L*v 21 # u = L*v
17 function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T 22 function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
18 uᵢ = L.a*apply_2nd_derivative(L.op, inverse_spacing(L.grid)[1], v, I[1]) 23 return apply(L.tensorOps[1],v,I)
19 return uᵢ
20 end 24 end
21 25
22 26
23 @inline function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T 27 @inline function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
24 # 2nd x-derivative 28 # 2nd x-derivative
25 @inbounds vx = view(v, :, Int(I[2])) 29 @inbounds vx = view(v, :, Int(I[2]))
26 @inbounds uᵢ = L.a*apply_2nd_derivative(L.op, inverse_spacing(L.grid)[1], vx , I[1]) 30 @inbounds uᵢ = apply(L.tensorOps[1], vx , (I[1],)) #Tuple conversion here is ugly. How to do it? Should we use indexing here?
31
27 # 2nd y-derivative 32 # 2nd y-derivative
28 @inbounds vy = view(v, Int(I[1]), :) 33 @inbounds vy = view(v, Int(I[1]), :)
29 @inbounds uᵢ += L.a*apply_2nd_derivative(L.op, inverse_spacing(L.grid)[2], vy, I[2]) 34 @inbounds uᵢ += apply(L.tensorOps[2], vy , (I[2],)) #Tuple conversion here is ugly. How to do it?
35
30 return uᵢ 36 return uᵢ
31 end 37 end
32 38
33 quadrature(L::Laplace) = Quadrature(L.op, L.grid) 39 quadrature(L::Laplace) = Quadrature(L.op, L.grid)
34 inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid) 40 inverse_quadrature(L::Laplace) = InverseQuadrature(L.op, L.grid)
35 boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId) 41 boundary_value(L::Laplace, bId::CartesianBoundary) = BoundaryValue(L.op, L.grid, bId)
36 normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId) 42 normal_derivative(L::Laplace, bId::CartesianBoundary) = NormalDerivative(L.op, L.grid, bId)
37 boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId) 43 boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId)
38 export quadrature 44 export quadrature
39 45
46 # At the moment the grid property is used all over. It could possibly be removed if we implement all the 1D operators as TensorMappings
40 """ 47 """
41 Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} 48 Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
42 49
43 Implements the quadrature operator `H` of Dim dimension as a TensorMapping 50 Implements the quadrature operator `H` of Dim dimension as a TensorMapping
44 """ 51 """