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

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
line wrap: on
line diff
--- a/SbpOperators/src/laplace/laplace.jl	Tue Jun 23 17:31:50 2020 +0200
+++ b/SbpOperators/src/laplace/laplace.jl	Tue Jun 23 17:32:54 2020 +0200
@@ -2,8 +2,8 @@
     Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
 
 Implements the Laplace operator `L` in Dim dimensions as a tensor operator
-The multi-dimensional tensor operator simply consists of a tuple of the 1D
-Laplace tensor operator as defined by ConstantLaplaceOperator.
+The multi-dimensional tensor operator consists of a tuple of 1D SecondDerivative
+tensor operators.
 """
 struct Laplace{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
     D2::NTuple(Dim,SecondDerivative{T,N,M,K})
@@ -17,20 +17,23 @@
     error("not implemented")
 end
 
+function LazyTensors.apply_transpose(L::Laplace{Dim,T}, v::AbstractArray{T,Dim}, I::NTuple{Dim,Index}) where {T,Dim} = LazyTensors.apply(L, v, I)
+
 # u = L*v
 function LazyTensors.apply(L::Laplace{1,T}, v::AbstractVector{T}, I::NTuple{1,Index}) where T
-    return apply(L.D2[1],v,I)
+    @inbounds u = apply(L.D2[1],v,I)
+    return u
 end
 
 
 @inline function LazyTensors.apply(L::Laplace{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T
     # 2nd x-derivative
     @inbounds vx = view(v, :, Int(I[2]))
-    @inbounds uᵢ = apply(L.D2[1], vx , (I[1],)) #Tuple conversion here is ugly. How to do it? Should we use indexing here?
+    @inbounds uᵢ = apply(L.D2[1], vx , I[1])
 
     # 2nd y-derivative
     @inbounds vy = view(v, Int(I[1]), :)
-    @inbounds uᵢ += apply(L.D2[2], vy , (I[2],)) #Tuple conversion here is ugly. How to do it?
+    @inbounds uᵢ += apply(L.D2[2], vy , I[2])
 
     return uᵢ
 end
@@ -42,31 +45,6 @@
 boundary_quadrature(L::Laplace, bId::CartesianBoundary) = BoundaryQuadrature(L.op, L.grid, bId)
 export quadrature
 
-# At the moment the grid property is used all over. It could possibly be removed if we implement all the 1D operators as TensorMappings
-"""
-    Quadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}
-
-Implements the quadrature operator `H` of Dim dimension as a TensorMapping
-"""
-struct Quadrature{Dim,T<:Real,N,M,K} <: TensorOperator{T,Dim}
-    op::D2{T,N,M,K}
-    grid::EquidistantGrid{Dim,T}
-end
-export Quadrature
-
-LazyTensors.domain_size(H::Quadrature{Dim}, range_size::NTuple{Dim,Integer}) where Dim = range_size
-
-@inline function LazyTensors.apply(H::Quadrature{2,T}, 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])
-    # Quadrature in y-direction
-    @inbounds q = apply_quadrature(H.op, spacing(H.grid)[2], q, I[2], N[2])
-    return q
-end
-
-LazyTensors.apply_transpose(H::Quadrature{2,T}, v::AbstractArray{T,2}, I::NTuple{2,Index}) where T = LazyTensors.apply(H,v,I)
-
 
 """
     InverseQuadrature{Dim,T<:Real,N,M,K} <: TensorMapping{T,Dim,Dim}