Mercurial > repos > public > sbplib_julia

comparison SbpOperators/src/laplace/laplace.jl @ 300:b00eea62c78e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create 1D tensor mapping for diagonal norm quadratures, and make the multi-dimensional quadrature use those. Move Qudrature from laplace.jl into Quadrature.jl
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 23 Jun 2020 17:32:54 +0200
parents f63232aeb1c6
children 6fa2ba769ae3
comparison
equal deleted inserted replaced
299:27a0bca5e1f2 300:b00eea62c78e
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 47
70 48
71 """ 49 """
72 InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim} 50 InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
73 51