text:
using CUDA
using Printf
using StaticArrays
using GPUifyLoops: @unroll

# include("deriv_ops.jl")

###
# y_in: input GPU array
# y_out: output GPU array
###

function D2x_GPU_naive(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM}) where {TILE_DIM}
	tidx = (blockIdx().x - 1) * TILE_DIM + threadIdx().x
	N = Nx*Ny
	# y_out = zeros(N)
	if tidx <= Ny
		y_out[tidx] = (y_in[tidx] - 2 * y_in[Ny + tidx] + y_in[2*Ny + tidx]) / h^2
	end
	sync_threads()

	if Ny+1 <= tidx <= N-Ny
		y_out[tidx] = (y_in[tidx - Ny] - 2 .* y_in[tidx] + y_in[tidx + Ny]) / h^2
	end

	sync_threads()

	if N-Ny+1 <= tidx <= N
		y_out[tidx] = (y_in[tidx - 2*Ny] -2 * y_in[tidx - Ny] + y_in[tidx]) / h^2
	end
	sync_threads()

	nothing
end

function D2x_GPU(y_in::CuArray,y_out::CuArray,Nx,Ny)
	y_out = similar(y_in)
	TILE_DIM_1 = 4
	TILE_DIM_2 = 16
	h = 1/Nx
	griddim = (div(Nx,TILE_DIM_1) + 1, div(Ny,TILE_DIM_2) + 1)
	blockdim = (TILE_DIM_1,TILE_DIM_2)
	@cuda threads=blockdim blocks=griddim D2x_GPU_shared(y_in, y_out, Nx, Ny, h, Val(TILE_DIM_1), Val(TILE_DIM_2))
	synchronize();
	return y_out
end

function D2x_GPU_shared(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
    tidx = threadIdx().x
    tidy = threadIdx().y

    # for global memory indexing
    i = (blockIdx().x - 1) * TILE_DIM1 + tidx
    j = (blockIdx().y - 1) * TILE_DIM2 + tidy

    global_index = i + (j - 1) * Ny

    HALO_WIDTH = 2 # For second order derivative

    tile = @cuStaticSharedMem(eltype(y_in), (TILE_DIM1, TILE_DIM2 + 2 * HALO_WIDTH))

    # for tile indexing
    k = tidx
    l = tidy

	# Writing pencil-shaped shared memory

	# for tile itself
	# if k <= TILE_DIM1 && l <= TILE_DIM2 && global_index <= Nx*Ny
	if k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny && j <= Nx
		# @inbounds tile[k,l+HALO_WIDTH] = y_in[global_index]
		 tile[k,l+HALO_WIDTH] = y_in[global_index]
    end

    sync_threads()

	# for left halo
	# if k <= TILE_DIM1 && l <= HALO_WIDTH && HALO_WIDTH*Ny+1 <= global_index <= (Nx+HALO_WIDTH)*Ny
	if k <= TILE_DIM1 && l <= HALO_WIDTH && i <= Ny && HALO_WIDTH+1 <= j <= HALO_WIDTH + Nx
		# @inbounds tile[k,l] = y_in[global_index - HALO_WIDTH*Ny]
		tile[k,l] = y_in[global_index - HALO_WIDTH*Ny]
	end

	sync_threads()

	# for right halo
	# if k <= TILE_DIM1 && l >= TILE_DIM2 - HALO_WIDTH && HALO_WIDTH*Ny+1 <= global_index <= (Nx-HALO_WIDTH)*Ny
	if k <= TILE_DIM1 && TILE_DIM2 - HALO_WIDTH + 1 <= l <= TILE_DIM2 && i <= Ny && j <= Nx - HALO_WIDTH
		# @inbounds tile[k,l+2*HALO_WIDTH] = y_in[global_index + HALO_WIDTH*Ny]
		tile[k,l+2*HALO_WIDTH] = y_in[global_index + HALO_WIDTH*Ny]
	end

    sync_threads()

    # Finite difference operation starts here

	# Left Boundary
	if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH -2 && i <= Ny && j == 1
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH] - 2*tile[k,l + HALO_WIDTH+1] + tile[k,l + HALO_WIDTH+2]) / h^2
		y_out[global_index] = (tile[k,l + HALO_WIDTH] - 2*tile[k,l + HALO_WIDTH+1] + tile[k,l + HALO_WIDTH+2]) / h^2
	end

	# Center
	if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH - 1 && i <= Ny && 2 <= j <= Nx-1
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-1] - 2*tile[k, l + HALO_WIDTH] + tile[k,l + HALO_WIDTH + 1]) / h^2
		y_out[global_index] = (tile[k,l + HALO_WIDTH-1] - 2*tile[k, l + HALO_WIDTH] + tile[k,l + HALO_WIDTH + 1]) / h^2
	end

	# Right Boundary
	if k <= TILE_DIM1 && 3 <= l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH && i <= Ny && j == Nx
		@inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-2] - 2*tile[k,l + HALO_WIDTH - 1] + tile[k,l + HALO_WIDTH]) / h^2
		# y_out[global_index] = (tile[k,l + HALO_WIDTH-2] - 2*tile[k,l + HALO_WIDTH - 1] + tile[k,l + HALO_WIDTH]) / h^2
		# y_out[global_index] = (tile[k,l+HALO_WIDTH-2])
		# y_out[global_index] = 0
	end

    sync_threads()

    nothing
end

function D2x_GPU_p4_naive(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM}) where {TILE_DIM}
	d = [1/12 -2/3 0 2/3 -1/12]
	bd = [-24/17  59/34  -4/17  -3/34  0     0;
	-1/2    0      1/2    0     0     0;
	 4/43 -59/86   0     59/86 -4/43  0;
	 3/98   0    -59/98   0    32/49 -4/49]
	tidx = (blockIdx().x - 1) * TILE_DIM + threadIdx().x
	N = Nx*Ny
	# y_out = zeros(N)
	if tidx <= Ny
		y_out[tidx] = (y_in[tidx] - 2 * y_in[Ny + tidx] + y_in[2*Ny + tidx]) / h^2
	end
	sync_threads()

	if Ny+1 <= tidx <= N-Ny
		y_out[tidx] = (y_in[tidx - Ny] - 2 .* y_in[tidx] + y_in[tidx + Ny]) / h^2
	end

	sync_threads()

	if N-Ny+1 <= tidx <= N
		y_out[tidx] = (y_in[tidx - 2*Ny] -2 * y_in[tidx - Ny] + y_in[tidx]) / h^2
	end
	sync_threads()

	nothing
end

function Dx_GPU(y_in::CuArray,y_out::CuArray,Nx,Ny)
	y_out = similar(y_in)
	TILE_DIM_1 = 4
	TILE_DIM_2 = 16
	h = 1/Nx
	griddim = (div(Nx,TILE_DIM_1) + 1, div(Ny,TILE_DIM_2) + 1)
	blockdim = (TILE_DIM_1,TILE_DIM_2)
	@cuda threads=blockdim blocks=griddim Dx_GPU_shared(y_in, y_out, Nx, Ny, h, Val(TILE_DIM_1), Val(TILE_DIM_2))
	synchronize();
	return y_out
end

function Dx_GPU_shared(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
	tidx = threadIdx().x
    tidy = threadIdx().y

    # for global memory indexing
    i = (blockIdx().x - 1) * TILE_DIM1 + tidx
    j = (blockIdx().y - 1) * TILE_DIM2 + tidy

    global_index = i + (j - 1) * Ny

    HALO_WIDTH = 2 # For second order derivative

    tile = @cuStaticSharedMem(eltype(y_in), (TILE_DIM1, TILE_DIM2 + 2 * HALO_WIDTH))

    # for tile indexing
    k = tidx
    l = tidy

	# Writing pencil-shaped shared memory

	# for tile itself
	# if k <= TILE_DIM1 && l <= TILE_DIM2 && global_index <= Nx*Ny
	if k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny && j <= Nx
		# @inbounds tile[k,l+HALO_WIDTH] = y_in[global_index]
		 tile[k,l+HALO_WIDTH] = y_in[global_index]
    end

    sync_threads()

	# for left halo
	# if k <= TILE_DIM1 && l <= HALO_WIDTH && HALO_WIDTH*Ny+1 <= global_index <= (Nx+HALO_WIDTH)*Ny
	if k <= TILE_DIM1 && l <= HALO_WIDTH && i <= Ny && HALO_WIDTH+1 <= j <= HALO_WIDTH + Nx
		# @inbounds tile[k,l] = y_in[global_index - HALO_WIDTH*Ny]
		tile[k,l] = y_in[global_index - HALO_WIDTH*Ny]
	end

	sync_threads()

	# for right halo
	# if k <= TILE_DIM1 && l >= TILE_DIM2 - HALO_WIDTH && HALO_WIDTH*Ny+1 <= global_index <= (Nx-HALO_WIDTH)*Ny
	if k <= TILE_DIM1 && TILE_DIM2 - HALO_WIDTH + 1<= l <= TILE_DIM2 && i <= Ny && j <= Nx - HALO_WIDTH
		# @inbounds tile[k,l+2*HALO_WIDTH] = y_in[global_index + HALO_WIDTH*Ny]
		tile[k,l+2*HALO_WIDTH] = y_in[global_index + HALO_WIDTH*Ny]
	end

    sync_threads()

    # Finite difference operation starts here

	# Left Boundary
	if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH -2 && i <= Ny && j == 1
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH] - 2*tile[k,l + HALO_WIDTH+1] + tile[k,l + HALO_WIDTH+2]) / h^2
		y_out[global_index] = (tile[k,l + HALO_WIDTH+1] - tile[k,l + HALO_WIDTH]) / h
	end

	# Center
	if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH - 1 && i <= Ny && 2 <= j <= Nx-1
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-1] - 2*tile[k, l + HALO_WIDTH] + tile[k,l + HALO_WIDTH + 1]) / h^2
		y_out[global_index] = (tile[k,l + HALO_WIDTH + 1] - tile[k,l + HALO_WIDTH - 1]) / (2*h)
	end

	# Right Boundary
	if k <= TILE_DIM1 && 3 <= l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH && i <= Ny && j == Nx
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-2] - 2*tile[k,l + HALO_WIDTH - 1] + tile[k,l + HALO_WIDTH]) / h^2
		y_out[global_index] = (tile[k,l + HALO_WIDTH ] - tile[k,l + HALO_WIDTH - 1]) / h
	end

    sync_threads()

    nothing
end

function Hxinv_GPU_shared(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
	tidx = threadIdx().x
    tidy = threadIdx().y

    # for global memory indexing
    i = (blockIdx().x - 1) * TILE_DIM1 + tidx
    j = (blockIdx().y - 1) * TILE_DIM2 + tidy

    global_index = i + (j - 1) * Ny

    HALO_WIDTH = 2 # For second order derivative

    tile = @cuStaticSharedMem(eltype(y_in), (TILE_DIM1, TILE_DIM2 + 2 * HALO_WIDTH))

    # for tile indexing
    k = tidx
    l = tidy

	# Writing pencil-shaped shared memory

	# for tile itself
	# if k <= TILE_DIM1 && l <= TILE_DIM2 && global_index <= Nx*Ny
	if k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny && j <= Nx
		# @inbounds tile[k,l+HALO_WIDTH] = y_in[global_index]
		 tile[k,l+HALO_WIDTH] = y_in[global_index]
    end

    sync_threads()

	# for left halo
	# if k <= TILE_DIM1 && l <= HALO_WIDTH && HALO_WIDTH*Ny+1 <= global_index <= (Nx+HALO_WIDTH)*Ny
	if k <= TILE_DIM1 && l <= HALO_WIDTH && i <= Ny && HALO_WIDTH+1 <= j <= HALO_WIDTH + Nx
		# @inbounds tile[k,l] = y_in[global_index - HALO_WIDTH*Ny]
		tile[k,l] = y_in[global_index - HALO_WIDTH*Ny]
	end

	sync_threads()

	# for right halo
	# if k <= TILE_DIM1 && l >= TILE_DIM2 - HALO_WIDTH && HALO_WIDTH*Ny+1 <= global_index <= (Nx-HALO_WIDTH)*Ny
	if k <= TILE_DIM1 && TILE_DIM2 - HALO_WIDTH + 1<= l <= TILE_DIM2 && i <= Ny && j <= Nx - HALO_WIDTH
		# @inbounds tile[k,l+2*HALO_WIDTH] = y_in[global_index + HALO_WIDTH*Ny]
		tile[k,l+2*HALO_WIDTH] = y_in[global_index + HALO_WIDTH*Ny]
	end

    sync_threads()

    # Finite difference operation starts here

	# Left Boundary
	if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH -2 && i <= Ny && j == 1
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH] - 2*tile[k,l + HALO_WIDTH+1] + tile[k,l + HALO_WIDTH+2]) / h^2
		y_out[global_index] = (2*tile[k,l + HALO_WIDTH]) / h
	end

	# Center
	if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH - 1 && i <= Ny && 2 <= j <= Nx-1
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-1] - 2*tile[k, l + HALO_WIDTH] + tile[k,l + HALO_WIDTH + 1]) / h^2
		y_out[global_index] = (tile[k,l + HALO_WIDTH] ) / (h)
	end

	# Right Boundary
	if k <= TILE_DIM1 && 3 <= l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH && i <= Ny && j == Nx
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-2] - 2*tile[k,l + HALO_WIDTH - 1] + tile[k,l + HALO_WIDTH]) / h^2
		y_out[global_index] = (2*tile[k,l + HALO_WIDTH ]) / h
	end

    sync_threads()

    nothing
end

function Hx_GPU_shared(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
	tidx = threadIdx().x
    tidy = threadIdx().y

    # for global memory indexing
    i = (blockIdx().x - 1) * TILE_DIM1 + tidx
    j = (blockIdx().y - 1) * TILE_DIM2 + tidy

    global_index = i + (j - 1) * Ny

    HALO_WIDTH = 2 # For second order derivative

    tile = @cuStaticSharedMem(eltype(y_in), (TILE_DIM1, TILE_DIM2 + 2 * HALO_WIDTH))

    # for tile indexing
    k = tidx
    l = tidy

	# Writing pencil-shaped shared memory

	# for tile itself
	# if k <= TILE_DIM1 && l <= TILE_DIM2 && global_index <= Nx*Ny
	if k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny && j <= Nx
		# @inbounds tile[k,l+HALO_WIDTH] = y_in[global_index]
		 tile[k,l+HALO_WIDTH] = y_in[global_index]
    end

    sync_threads()

	# for left halo
	# if k <= TILE_DIM1 && l <= HALO_WIDTH && HALO_WIDTH*Ny+1 <= global_index <= (Nx+HALO_WIDTH)*Ny
	if k <= TILE_DIM1 && l <= HALO_WIDTH && i <= Ny && HALO_WIDTH+1 <= j <= HALO_WIDTH + Nx
		# @inbounds tile[k,l] = y_in[global_index - HALO_WIDTH*Ny]
		tile[k,l] = y_in[global_index - HALO_WIDTH*Ny]
	end

	sync_threads()

	# for right halo
	# if k <= TILE_DIM1 && l >= TILE_DIM2 - HALO_WIDTH && HALO_WIDTH*Ny+1 <= global_index <= (Nx-HALO_WIDTH)*Ny
	if k <= TILE_DIM1 && TILE_DIM2 - HALO_WIDTH + 1 <= l <= TILE_DIM2 && i <= Ny && j <= Nx - HALO_WIDTH
		# @inbounds tile[k,l+2*HALO_WIDTH] = y_in[global_index + HALO_WIDTH*Ny]
		tile[k,l+2*HALO_WIDTH] = y_in[global_index + HALO_WIDTH*Ny]
	end

    sync_threads()

    # Finite difference operation starts here

	# Left Boundary
	if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH -2 && i <= Ny && j == 1
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH] - 2*tile[k,l + HALO_WIDTH+1] + tile[k,l + HALO_WIDTH+2]) / h^2
		y_out[global_index] = (h*tile[k,l + HALO_WIDTH]) / 2
	end

	# Center
	if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH - 1 && i <= Ny && 2 <= j <= Nx-1
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-1] - 2*tile[k, l + HALO_WIDTH] + tile[k,l + HALO_WIDTH + 1]) / h^2
		y_out[global_index] = h*(tile[k,l + HALO_WIDTH] )
	end

	# Right Boundary
	if k <= TILE_DIM1 && 3 <= l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH && i <= Ny && j == Nx
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-2] - 2*tile[k,l + HALO_WIDTH - 1] + tile[k,l + HALO_WIDTH]) / h^2
		y_out[global_index] = (h*tile[k,l + HALO_WIDTH ]) / 2
	end

    sync_threads()

    nothing
end

function D2y_GPU_naive(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM}) where {TILE_DIM}
	tidx = (blockIdx().x - 1) * TILE_DIM + threadIdx().x
	N = Nx*Ny
	if 2 <= tidx <= N-1
		@inbounds y_out[tidx] = (y_in[tidx-1] - 2y_in[tidx] + y_in[tidx + 1]) / h^2
	end

	if 1 <= tidx <= N && mod(tidx,Ny) == 0
		@inbounds y_out[tidx] = (y_in[tidx] - 2y_in[tidx - 1] + y_in[tidx - 2]) / h^2
		@inbounds y_out[tidx-Ny+1] = (y_in[tidx-Ny+1] - 2y_in[tidx - Ny + 2] + y_in[tidx - Ny + 3]) / h^2
	end

	sync_threads()
	nothing
end

function D2y_GPU_shared(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
    tidx = threadIdx().x
    tidy = threadIdx().y

	i = (blockIdx().x - 1) * TILE_DIM1 + tidx
	j = (blockIdx().y - 1) * TILE_DIM2 + tidy

	global_index = i + (j-1)*Nx

	HALO_WIDTH = 2
	tile = @cuStaticSharedMem(eltype(y_in),(TILE_DIM1+2*HALO_WIDTH,TILE_DIM2))

	k = tidx
	l = tidy

    # Writing pencil-shaped shared memory

    # for tile itself
	# if k <= TILE_DIM1 && l <= TILE_DIM2 && global_index <= Nx*Ny
	if k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny && j <= Nx
		# @inbounds tile[k+HALO_WIDTH,l] = y_in[global_index]
		tile[k+HALO_WIDTH,l] = y_in[global_index]
	end

	sync_threads()

	# For upper halo
	# if k <= HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= global_index <= Nx*Ny + HALO_WIDTH
	if k <= HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= i <= Ny + HALO_WIDTH && j <= Nx
		# @inbounds tile[k,l] = y_in[global_index - HALO_WIDTH]
		tile[k,l] = y_in[global_index - HALO_WIDTH]
	end

	sync_threads()

	# For lower halo
	# if k >= TILE_DIM1 - HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= global_index <= Nx*Ny - HALO_WIDTH
	if  TILE_DIM1 - HALO_WIDTH + 1 <= k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny - HALO_WIDTH && j <= Nx
		# @inbounds tile[k+2*HALO_WIDTH,l] = y_in[global_index + HALO_WIDTH]
		tile[k+2*HALO_WIDTH,l] = y_in[global_index + HALO_WIDTH]
	end

    sync_threads()

    # Finite Difference Operations starts

    #Upper Boundary
	if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH -2 && l <= TILE_DIM2 && i == 1 && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH,l] - 2*tile[k+HALO_WIDTH+1,l] + tile[k+HALO_WIDTH+2,l]) / h^2
		y_out[global_index] = (tile[k+HALO_WIDTH,l] - 2*tile[k+HALO_WIDTH+1,l] + tile[k+HALO_WIDTH+2,l]) / h^2
	end

	sync_threads()

	#Center
	if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH - 1 && l <= TILE_DIM2 && 2 <= i <= Nx-1 && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH-1,l] - 2*tile[k+HALO_WIDTH,l] + tile[k+HALO_WIDTH+1,l]) / h^2
		y_out[global_index] = (tile[k+HALO_WIDTH-1,l] - 2*tile[k+HALO_WIDTH,l] + tile[k+HALO_WIDTH+1,l]) / h^2
	end

	sync_threads()

	#Lower Boundary
	if 3 <= k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && i == Nx && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH-2,l] - 2*tile[k+HALO_WIDTH-1,l] + tile[k+HALO_WIDTH,l]) / h^2
		y_out[global_index] = (tile[k+HALO_WIDTH-2,l] - 2*tile[k+HALO_WIDTH-1,l] + tile[k+HALO_WIDTH,l]) / h^2
		# y_out[global_index] = (tile[k+HALO_WIDTH-2,l] - 2*tile[k+HALO_WIDTH-1,l])
		# y_out[global_index] = tile[k+HALO_WIDTH,l]
		# y_out[global_index] = 0
	end

	# if i == 1 && j <= Nx
	# 	y_out[global_index] = y_out[global_index+1]
	# end

	# sync_threads()

	# if i == Ny && j <= Nx
	# 	y_out[global_index] = y_out[global_index-1]
	# 	# y_out[global_index] = 0.0
	# end

    sync_threads()

    nothing

end

function Dy_GPU_shared(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
	tidx = threadIdx().x
    tidy = threadIdx().y

	i = (blockIdx().x - 1) * TILE_DIM1 + tidx
	j = (blockIdx().y - 1) * TILE_DIM2 + tidy

	global_index = i + (j-1)*Nx

	HALO_WIDTH = 2
	tile = @cuStaticSharedMem(eltype(y_in),(TILE_DIM1+2*HALO_WIDTH,TILE_DIM2))

	k = tidx
	l = tidy

    # Writing pencil-shaped shared memory

    # for tile itself
	# if k <= TILE_DIM1 && l <= TILE_DIM2 && global_index <= Nx*Ny
	if k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny && j <= Nx
		# @inbounds tile[k+HALO_WIDTH,l] = y_in[global_index]
		tile[k+HALO_WIDTH,l] = y_in[global_index]
	end

	sync_threads()

	# For upper halo
	# if k <= HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= global_index <= Nx*Ny + HALO_WIDTH
	if k <= HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= i <= Ny + HALO_WIDTH && j <= Nx
		# @inbounds tile[k,l] = y_in[global_index - HALO_WIDTH]
		tile[k,l] = y_in[global_index - HALO_WIDTH]
	end

	sync_threads()

	# For lower halo
	# if k >= TILE_DIM1 - HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= global_index <= Nx*Ny - HALO_WIDTH
	if  TILE_DIM1 - HALO_WIDTH + 1 <= k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny - HALO_WIDTH && j <= Nx
		# @inbounds tile[k+2*HALO_WIDTH,l] = y_in[global_index + HALO_WIDTH]
		tile[k+2*HALO_WIDTH,l] = y_in[global_index + HALO_WIDTH]
	end

    sync_threads()

    # Finite Difference Operations starts

    #Upper Boundary
	if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH - 2 && l <= TILE_DIM2 && i == 1 && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH,l] - 2*tile[k+HALO_WIDTH+1,l] + tile[k+HALO_WIDTH+2,l]) / h^2
		y_out[global_index] = (tile[k+HALO_WIDTH+1,l] - tile[k+HALO_WIDTH,l]) / h
	end

	sync_threads()

	#Center
	if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && 2 <= i <= Nx-1 && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH-1,l] - 2*tile[k+HALO_WIDTH,l] + tile[k+HALO_WIDTH+1,l]) / h^2
		y_out[global_index] = (tile[k+HALO_WIDTH+1,l] - tile[k+HALO_WIDTH-1,l]) / (2*h)
	end

	sync_threads()

	#Lower Boundary
	if 3 <= k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && i == Nx && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH-2,l] - 2*tile[k+HALO_WIDTH-1,l] + tile[k+HALO_WIDTH,l]) / h^2
		y_out[global_index] = (tile[k+HALO_WIDTH,l] - tile[k+HALO_WIDTH-1,l]) / h
		# y_out[global_index] = tile[k+HALO_WIDTH-1,l] #- tile[k+HALO_WIDTH-1,l]
 	end

    sync_threads()

    nothing
end

function Hyinv_GPU_shared(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
	tidx = threadIdx().x
    tidy = threadIdx().y

	i = (blockIdx().x - 1) * TILE_DIM1 + tidx
	j = (blockIdx().y - 1) * TILE_DIM2 + tidy

	global_index = i + (j-1)*Nx

	HALO_WIDTH = 2
	tile = @cuStaticSharedMem(eltype(y_in),(TILE_DIM1+2*HALO_WIDTH,TILE_DIM2))

	k = tidx
	l = tidy

    # Writing pencil-shaped shared memory

    # for tile itself
	# if k <= TILE_DIM1 && l <= TILE_DIM2 && global_index <= Nx*Ny
	if k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny && j <= Nx
		# @inbounds tile[k+HALO_WIDTH,l] = y_in[global_index]
		tile[k+HALO_WIDTH,l] = y_in[global_index]
	end

	sync_threads()

	# For upper halo
	# if k <= HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= global_index <= Nx*Ny + HALO_WIDTH
	if k <= HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= i <= Ny + HALO_WIDTH && j <= Nx
		# @inbounds tile[k,l] = y_in[global_index - HALO_WIDTH]
		tile[k,l] = y_in[global_index - HALO_WIDTH]
	end

	sync_threads()

	# For lower halo
	# if k >= TILE_DIM1 - HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= global_index <= Nx*Ny - HALO_WIDTH
	if  TILE_DIM1 - HALO_WIDTH <= k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny - HALO_WIDTH && j <= Nx
		# @inbounds tile[k+2*HALO_WIDTH,l] = y_in[global_index + HALO_WIDTH]
		tile[k+2*HALO_WIDTH,l] = y_in[global_index + HALO_WIDTH]
	end

    sync_threads()

    # Finite Difference Operations starts

    #Upper Boundary
	if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && i == 1 && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH,l] - 2*tile[k+HALO_WIDTH+1,l] + tile[k+HALO_WIDTH+2,l]) / h^2
		y_out[global_index] = (2*tile[k+HALO_WIDTH,l]) / h
	end

	sync_threads()

	#Center
	if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && 2 <= i <= Nx-1 && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH-1,l] - 2*tile[k+HALO_WIDTH,l] + tile[k+HALO_WIDTH+1,l]) / h^2
		y_out[global_index] = (tile[k+HALO_WIDTH,l]) / h
	end

	sync_threads()

	#Lower Boundary
	if 3 <= k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && i == Nx && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH-2,l] - 2*tile[k+HALO_WIDTH-1,l] + tile[k+HALO_WIDTH,l]) / h^2
		y_out[global_index] = (2*tile[k+HALO_WIDTH,l]) / h
    end

    sync_threads()

    nothing
end

function Hy_GPU_shared(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
	tidx = threadIdx().x
    tidy = threadIdx().y

	i = (blockIdx().x - 1) * TILE_DIM1 + tidx
	j = (blockIdx().y - 1) * TILE_DIM2 + tidy

	global_index = i + (j-1)*Ny

	HALO_WIDTH = 2
	tile = @cuStaticSharedMem(eltype(y_in),(TILE_DIM1+2*HALO_WIDTH,TILE_DIM2))

	k = tidx
	l = tidy

    # Writing pencil-shaped shared memory

    # for tile itself
	# if k <= TILE_DIM1 && l <= TILE_DIM2 && global_index <= Nx*Ny
	if k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny && j <= Nx
		# @inbounds tile[k+HALO_WIDTH,l] = y_in[global_index]
		tile[k+HALO_WIDTH,l] = y_in[global_index]
	end

	sync_threads()

	# For upper halo
	# if k <= HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= global_index <= Nx*Ny + HALO_WIDTH
	if k <= HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= i <= Ny + HALO_WIDTH && j <= Nx
		# @inbounds tile[k,l] = y_in[global_index - HALO_WIDTH]
		tile[k,l] = y_in[global_index - HALO_WIDTH]
	end

	sync_threads()

	# For lower halo
	# if k >= TILE_DIM1 - HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= global_index <= Nx*Ny - HALO_WIDTH
	if  TILE_DIM1 - HALO_WIDTH + 1 <= k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny - HALO_WIDTH && j <= Nx
		# @inbounds tile[k+2*HALO_WIDTH,l] = y_in[global_index + HALO_WIDTH]
		tile[k+2*HALO_WIDTH,l] = y_in[global_index + HALO_WIDTH]
	end

    sync_threads()

    # Finite Difference Operations starts

    #Upper Boundary
	# if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH -2 && l <= TILE_DIM2 && i == 1 && j <= Ny
	if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && i == 1 && j <= Nx
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH,l] - 2*tile[k+HALO_WIDTH+1,l] + tile[k+HALO_WIDTH+2,l]) / h^2
		y_out[global_index] = (h*tile[k+HALO_WIDTH,l]) / 2
	end

	# sync_threads()

	#Center
	# if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH - 1 && l <= TILE_DIM2 && 2 <= i <= Nx-1 && j <= Ny
	if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && 2 <= i <= Ny-1 && j <= Nx
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH-1,l] - 2*tile[k+HALO_WIDTH,l] + tile[k+HALO_WIDTH+1,l]) / h^2
		y_out[global_index] = h * (tile[k+HALO_WIDTH,l])
	end

	# sync_threads()

	#Lower Boundary
	if 1 <= k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && i == Ny && j <= Nx
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH-2,l] - 2*tile[k+HALO_WIDTH-1,l] + tile[k+HALO_WIDTH,l]) / h^2
		y_out[global_index] = (h*tile[k+HALO_WIDTH,l]) / 2
    end

    sync_threads()

    nothing
end

function Bx_GPU_shared(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
	tidx = threadIdx().x
    tidy = threadIdx().y

    # for global memory indexing
    i = (blockIdx().x - 1) * TILE_DIM1 + tidx
    j = (blockIdx().y - 1) * TILE_DIM2 + tidy

    global_index = i + (j - 1) * Ny

    HALO_WIDTH = 2 # For second order derivative

    tile = @cuStaticSharedMem(eltype(y_in), (TILE_DIM1, TILE_DIM2 + 2 * HALO_WIDTH))

    # for tile indexing
    k = tidx
    l = tidy

	# Writing pencil-shaped shared memory

	# for tile itself
	# if k <= TILE_DIM1 && l <= TILE_DIM2 && global_index <= Nx*Ny
	if k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny && j <= Nx
		# @inbounds tile[k,l+HALO_WIDTH] = y_in[global_index]
		 tile[k,l+HALO_WIDTH] = y_in[global_index]
    end

    sync_threads()

	# for left halo
	# if k <= TILE_DIM1 && l <= HALO_WIDTH && HALO_WIDTH*Ny+1 <= global_index <= (Nx+HALO_WIDTH)*Ny
	if k <= TILE_DIM1 && l <= HALO_WIDTH && i <= Ny && HALO_WIDTH+1 <= j <= HALO_WIDTH + Nx
		# @inbounds tile[k,l] = y_in[global_index - HALO_WIDTH*Ny]
		tile[k,l] = y_in[global_index - HALO_WIDTH*Ny]
	end

	sync_threads()

	# for right halo
	# if k <= TILE_DIM1 && l >= TILE_DIM2 - HALO_WIDTH && HALO_WIDTH*Ny+1 <= global_index <= (Nx-HALO_WIDTH)*Ny
	if k <= TILE_DIM1 && TILE_DIM2 - HALO_WIDTH + 1<= l <= TILE_DIM2 && i <= Ny && j <= Nx - HALO_WIDTH
		# @inbounds tile[k,l+2*HALO_WIDTH] = y_in[global_index + HALO_WIDTH*Ny]
		tile[k,l+2*HALO_WIDTH] = y_in[global_index + HALO_WIDTH*Ny]
	end

    sync_threads()

    # Finite difference operation starts here

	# Left Boundary
	if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH -2 && i <= Ny && j == 1
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH] - 2*tile[k,l + HALO_WIDTH+1] + tile[k,l + HALO_WIDTH+2]) / h^2
		y_out[global_index] = -1.0
	end

	# Center
	if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH - 1 && i <= Ny && 2 <= j <= Nx-1
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-1] - 2*tile[k, l + HALO_WIDTH] + tile[k,l + HALO_WIDTH + 1]) / h^2
		y_out[global_index] = 0.0
	end

	# Right Boundary
	if k <= TILE_DIM1 && 3 <= l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH && i <= Ny && j == Nx
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-2] - 2*tile[k,l + HALO_WIDTH - 1] + tile[k,l + HALO_WIDTH]) / h^2
		y_out[global_index] = 1.0
	end

    sync_threads()

    nothing
end

function BxSx_GPU_shared(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
	tidx = threadIdx().x
    tidy = threadIdx().y

    # for global memory indexing
    i = (blockIdx().x - 1) * TILE_DIM1 + tidx
    j = (blockIdx().y - 1) * TILE_DIM2 + tidy

    global_index = i + (j - 1) * Ny

    HALO_WIDTH = 2 # For second order derivative

    tile = @cuStaticSharedMem(eltype(y_in), (TILE_DIM1, TILE_DIM2 + 2 * HALO_WIDTH))

    # for tile indexing
    k = tidx
    l = tidy

	# Writing pencil-shaped shared memory

	# for tile itself
	# if k <= TILE_DIM1 && l <= TILE_DIM2 && global_index <= Nx*Ny
	if k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny && j <= Nx
		# @inbounds tile[k,l+HALO_WIDTH] = y_in[global_index]
		 tile[k,l+HALO_WIDTH] = y_in[global_index]
    end

    sync_threads()

	# for left halo
	# if k <= TILE_DIM1 && l <= HALO_WIDTH && HALO_WIDTH*Ny+1 <= global_index <= (Nx+HALO_WIDTH)*Ny
	if k <= TILE_DIM1 && l <= HALO_WIDTH && i <= Ny && HALO_WIDTH+1 <= j <= HALO_WIDTH + Nx
		# @inbounds tile[k,l] = y_in[global_index - HALO_WIDTH*Ny]
		tile[k,l] = y_in[global_index - HALO_WIDTH*Ny]
	end

	sync_threads()

	# for right halo
	# if k <= TILE_DIM1 && l >= TILE_DIM2 - HALO_WIDTH && HALO_WIDTH*Ny+1 <= global_index <= (Nx-HALO_WIDTH)*Ny
	if k <= TILE_DIM1 && TILE_DIM2 - HALO_WIDTH + 1 <= l <= TILE_DIM2 && i <= Ny && j <= Nx - HALO_WIDTH
		# @inbounds tile[k,l+2*HALO_WIDTH] = y_in[global_index + HALO_WIDTH*Ny]
		tile[k,l+2*HALO_WIDTH] = y_in[global_index + HALO_WIDTH*Ny]
	end

    sync_threads()

    # Finite difference operation starts here

	# Left Boundary
	if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH -2 && i <= Ny && j == 1
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH] - 2*tile[k,l + HALO_WIDTH+1] + tile[k,l + HALO_WIDTH+2]) / h^2
		y_out[global_index] = (1.5 * tile[k,l+HALO_WIDTH] - 2.0 * tile[k,l+HALO_WIDTH+1] + 0.5*tile[k,l+HALO_WIDTH+2]) / h
	end

	# Center
	if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH - 1 && i <= Ny && 2 <= j <= Nx-1
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-1] - 2*tile[k, l + HALO_WIDTH] + tile[k,l + HALO_WIDTH + 1]) / h^2
		y_out[global_index] = 0.0
	end

	# Right Boundary
	if k <= TILE_DIM1 && 3 <= l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH && i <= Ny && j == Nx
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-2] - 2*tile[k,l + HALO_WIDTH - 1] + tile[k,l + HALO_WIDTH]) / h^2
		y_out[global_index] = (0.5 * tile[k,l+HALO_WIDTH-2] - 2.0 * tile[k,l+HALO_WIDTH-1] + 1.5 * tile[k,l+HALO_WIDTH]) / h
	end

    sync_threads()

    nothing
end

function By_GPU_shared(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
	tidx = threadIdx().x
    tidy = threadIdx().y

	i = (blockIdx().x - 1) * TILE_DIM1 + tidx
	j = (blockIdx().y - 1) * TILE_DIM2 + tidy

	global_index = i + (j-1)*Nx

	HALO_WIDTH = 2
	tile = @cuStaticSharedMem(eltype(y_in),(TILE_DIM1+2*HALO_WIDTH,TILE_DIM2))

	k = tidx
	l = tidy

    # Writing pencil-shaped shared memory

    # for tile itself
	# if k <= TILE_DIM1 && l <= TILE_DIM2 && global_index <= Nx*Ny
	if k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny && j <= Nx
		# @inbounds tile[k+HALO_WIDTH,l] = y_in[global_index]
		tile[k+HALO_WIDTH,l] = y_in[global_index]
	end

	sync_threads()

	# For upper halo
	# if k <= HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= global_index <= Nx*Ny + HALO_WIDTH
	if k <= HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= i <= Ny && j <= Nx
		# @inbounds tile[k,l] = y_in[global_index - HALO_WIDTH]
		tile[k,l] = y_in[global_index - HALO_WIDTH]
	end

	sync_threads()

	# For lower halo
	# if k >= TILE_DIM1 - HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= global_index <= Nx*Ny - HALO_WIDTH
	if  TILE_DIM1 - HALO_WIDTH <= k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny - HALO_WIDTH && j <= Nx
		# @inbounds tile[k+2*HALO_WIDTH,l] = y_in[global_index + HALO_WIDTH]
		tile[k+2*HALO_WIDTH,l] = y_in[global_index + HALO_WIDTH]
	end

    sync_threads()

    # Finite Difference Operations starts

    #Upper Boundary
	if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && i == 1 && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH,l] - 2*tile[k+HALO_WIDTH+1,l] + tile[k+HALO_WIDTH+2,l]) / h^2
		y_out[global_index] = -1.0
	end

	sync_threads()

	#Center
	if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && 2 <= i <= Nx-1 && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH-1,l] - 2*tile[k+HALO_WIDTH,l] + tile[k+HALO_WIDTH+1,l]) / h^2
		y_out[global_index] = 0.0
	end

	sync_threads()

	#Lower Boundary
	if 3 <= k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && i == Nx && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH-2,l] - 2*tile[k+HALO_WIDTH-1,l] + tile[k+HALO_WIDTH,l]) / h^2
		y_out[global_index] = 1.0
    end

    sync_threads()

	nothing
end

function BySy_GPU_shared(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
    tidx = threadIdx().x
    tidy = threadIdx().y

	i = (blockIdx().x - 1) * TILE_DIM1 + tidx
	j = (blockIdx().y - 1) * TILE_DIM2 + tidy

	global_index = i + (j-1)*Nx

	HALO_WIDTH = 2
	tile = @cuStaticSharedMem(eltype(y_in),(TILE_DIM1+2*HALO_WIDTH,TILE_DIM2))

	k = tidx
	l = tidy

    # Writing pencil-shaped shared memory

    # for tile itself
	# if k <= TILE_DIM1 && l <= TILE_DIM2 && global_index <= Nx*Ny
	if k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny && j <= Nx
		# @inbounds tile[k+HALO_WIDTH,l] = y_in[global_index]
		tile[k+HALO_WIDTH,l] = y_in[global_index]
	end

	sync_threads()

	# For upper halo
	# if k <= HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= global_index <= Nx*Ny + HALO_WIDTH
	if k <= HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= i <= Ny + HALO_WIDTH && j <= Nx
		# @inbounds tile[k,l] = y_in[global_index - HALO_WIDTH]
		tile[k,l] = y_in[global_index - HALO_WIDTH]
	end

	sync_threads()

	# For lower halo
	# if k >= TILE_DIM1 - HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= global_index <= Nx*Ny - HALO_WIDTH
	if  TILE_DIM1 - HALO_WIDTH <= k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny - HALO_WIDTH && j <= Nx
		# @inbounds tile[k+2*HALO_WIDTH,l] = y_in[global_index + HALO_WIDTH]
		tile[k+2*HALO_WIDTH,l] = y_in[global_index + HALO_WIDTH]
	end

    sync_threads()

    # Finite Difference Operations starts

    #Upper Boundary
	if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && i == 1 && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH,l] - 2*tile[k+HALO_WIDTH+1,l] + tile[k+HALO_WIDTH+2,l]) / h^2
		y_out[global_index] = (1.5*tile[k+HALO_WIDTH,l] - 2*tile[k+HALO_WIDTH+1,l] + 0.5*tile[k+HALO_WIDTH+2,l]) / h
	end

	sync_threads()

	#Center
	if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && 2 <= i <= Nx-1 && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH-1,l] - 2*tile[k+HALO_WIDTH,l] + tile[k+HALO_WIDTH+1,l]) / h^2
		y_out[global_index] = 0.0
	end

	sync_threads()

	#Lower Boundary
	if 3 <= k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && i == Nx && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH-2,l] - 2*tile[k+HALO_WIDTH-1,l] + tile[k+HALO_WIDTH,l]) / h^2
		y_out[global_index] = (0.5 * tile[k+HALO_WIDTH-2,l] - 2*tile[k+HALO_WIDTH-1,l] + 1.5*tile[k+HALO_WIDTH,l]) / h
    end

    sync_threads()

    nothing

end

function BxSx_tran_GPU_shared(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
	tidx = threadIdx().x
    tidy = threadIdx().y

    # for global memory indexing
    i = (blockIdx().x - 1) * TILE_DIM1 + tidx
    j = (blockIdx().y - 1) * TILE_DIM2 + tidy

    global_index = i + (j - 1) * Ny

    HALO_WIDTH = 2 # For second order derivative

    tile = @cuStaticSharedMem(eltype(y_in), (TILE_DIM1, TILE_DIM2 + 2 * HALO_WIDTH))

    # for tile indexing
    k = tidx
    l = tidy

	# Writing pencil-shaped shared memory

	# for tile itself
	# if k <= TILE_DIM1 && l <= TILE_DIM2 && global_index <= Nx*Ny
	if k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny && j <= Nx
		# @inbounds tile[k,l+HALO_WIDTH] = y_in[global_index]
		 tile[k,l+HALO_WIDTH] = y_in[global_index]
    end

    sync_threads()

	# for left halo
	# if k <= TILE_DIM1 && l <= HALO_WIDTH && HALO_WIDTH*Ny+1 <= global_index <= (Nx+HALO_WIDTH)*Ny
	if k <= TILE_DIM1 && l <= HALO_WIDTH && i <= Ny && HALO_WIDTH+1 <= j <= HALO_WIDTH + Nx
		# @inbounds tile[k,l] = y_in[global_index - HALO_WIDTH*Ny]
		tile[k,l] = y_in[global_index - HALO_WIDTH*Ny]
	end

	sync_threads()

	# for right halo
	# if k <= TILE_DIM1 && l >= TILE_DIM2 - HALO_WIDTH && HALO_WIDTH*Ny+1 <= global_index <= (Nx-HALO_WIDTH)*Ny
	if k <= TILE_DIM1 && TILE_DIM2 - HALO_WIDTH <= l <= TILE_DIM2 && i <= Ny && j <= Nx - HALO_WIDTH
		# @inbounds tile[k,l+2*HALO_WIDTH] = y_in[global_index + HALO_WIDTH*Ny]
		tile[k,l+2*HALO_WIDTH] = y_in[global_index + HALO_WIDTH*Ny]
	end

    sync_threads()

	# Finite difference operation starts here

	# Center
	if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH - 1 && i <= Ny && 4 <= j <= Nx-3
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-1] - 2*tile[k, l + HALO_WIDTH] + tile[k,l + HALO_WIDTH + 1]) / h^2
		y_out[global_index] = 0.0
	end

	sync_threads()

	# Left Boundary
	if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH -2 && i <= Ny && j == 1
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH] - 2*tile[k,l + HALO_WIDTH+1] + tile[k,l + HALO_WIDTH+2]) / h^2
		y_out[global_index] = (1.5 * tile[k,l+HALO_WIDTH]) / h
		y_out[global_index + Ny] = (-2.0 * tile[k,l+HALO_WIDTH])/h
		y_out[global_index + 2*Ny] = (0.5 * tile[k,l+HALO_WIDTH])/h
	end
	sync_threads()

	# if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH -2 && i <= Ny && j == 2
	# 	# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH] - 2*tile[k,l + HALO_WIDTH+1] + tile[k,l + HALO_WIDTH+2]) / h^2
	# 	y_out[global_index] = (-2.0 * tile[k,l+HALO_WIDTH]) / h
	# end

	# if k <= TILE_DIM1 && l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH -2 && i <= Ny && j == 3
	# 	# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH] - 2*tile[k,l + HALO_WIDTH+1] + tile[k,l + HALO_WIDTH+2]) / h^2
	# 	y_out[global_index] = (0.5 * tile[k,l+HALO_WIDTH]) / h
	# end

	# Right Boundary
	if k <= TILE_DIM1 && 3 <= l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH && i <= Ny && j == Nx
		# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-2] - 2*tile[k,l + HALO_WIDTH - 1] + tile[k,l + HALO_WIDTH]) / h^2
		y_out[global_index] = (1.5 * tile[k,l+HALO_WIDTH]) / h
		y_out[global_index - Ny] = (-2.0 * tile[k, l+HALO_WIDTH])/h
		y_out[global_index - 2*Ny] = (0.5 * tile[k, l+HALO_WIDTH])/h
	end
	sync_threads()

	# if k <= TILE_DIM1 && 3 <= l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH && i <= Ny && j == Nx - 1
	# 	# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-2] - 2*tile[k,l + HALO_WIDTH - 1] + tile[k,l + HALO_WIDTH]) / h^2
	# 	y_out[global_index] = (-2.0 * tile[k,l+HALO_WIDTH]) / h
	# end

	# if k <= TILE_DIM1 && 3 <= l + HALO_WIDTH <= TILE_DIM2 + 2*HALO_WIDTH && i <= Ny && j == Nx - 2
	# 	# @inbounds y_out[global_index] = (tile[k,l + HALO_WIDTH-2] - 2*tile[k,l + HALO_WIDTH - 1] + tile[k,l + HALO_WIDTH]) / h^2
	# 	y_out[global_index] = (0.5 * tile[k,l+HALO_WIDTH]) / h
	# end

    nothing
end

function BySy_tran_GPU_shared(y_in, y_out, Nx, Ny, h, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
    tidx = threadIdx().x
    tidy = threadIdx().y

	i = (blockIdx().x - 1) * TILE_DIM1 + tidx
	j = (blockIdx().y - 1) * TILE_DIM2 + tidy

	global_index = i + (j-1)*Nx

	HALO_WIDTH = 2
	tile = @cuStaticSharedMem(eltype(y_in),(TILE_DIM1+2*HALO_WIDTH,TILE_DIM2))

	k = tidx
	l = tidy

    # Writing pencil-shaped shared memory

    # for tile itself
	# if k <= TILE_DIM1 && l <= TILE_DIM2 && global_index <= Nx*Ny
	if k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny && j <= Nx
		# @inbounds tile[k+HALO_WIDTH,l] = y_in[global_index]
		tile[k+HALO_WIDTH,l] = y_in[global_index]
	end

	sync_threads()

	# For upper halo
	# if k <= HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= global_index <= Nx*Ny + HALO_WIDTH
	if k <= HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= i <= Ny && j <= Nx
		# @inbounds tile[k,l] = y_in[global_index - HALO_WIDTH]
		tile[k,l] = y_in[global_index - HALO_WIDTH]
	end

	sync_threads()

	# For lower halo
	# if k >= TILE_DIM1 - HALO_WIDTH && l <= TILE_DIM2 && HALO_WIDTH + 1 <= global_index <= Nx*Ny - HALO_WIDTH
	if  TILE_DIM1 - HALO_WIDTH <= k <= TILE_DIM1 && l <= TILE_DIM2 && i <= Ny - HALO_WIDTH && j <= Nx
		# @inbounds tile[k+2*HALO_WIDTH,l] = y_in[global_index + HALO_WIDTH]
		tile[k+2*HALO_WIDTH,l] = y_in[global_index + HALO_WIDTH]
	end

    sync_threads()

    # Finite Difference Operations starts

    #Upper Boundary
	if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && i == 1 && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH,l] - 2*tile[k+HALO_WIDTH+1,l] + tile[k+HALO_WIDTH+2,l]) / h^2
		y_out[global_index] = (1.5*tile[k+HALO_WIDTH,l]) / h
		y_out[global_index + 1] = (-2.0*tile[k+HALO_WIDTH,l]) / h
		y_out[global_index + 2] = (0.5*tile[k+HALO_WIDTH,l])/h
	end

	sync_threads()

	#Center
	if k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && 4 <= i <= Nx-3 && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH-1,l] - 2*tile[k+HALO_WIDTH,l] + tile[k+HALO_WIDTH+1,l]) / h^2
		y_out[global_index] = 0.0
	end

	sync_threads()

	#Lower Boundary
	if 3 <= k + HALO_WIDTH <= TILE_DIM1 + 2*HALO_WIDTH && l <= TILE_DIM2 && i == Nx && j <= Ny
		# @inbounds y_out[global_index] = (tile[k+HALO_WIDTH-2,l] - 2*tile[k+HALO_WIDTH-1,l] + tile[k+HALO_WIDTH,l]) / h^2
		y_out[global_index] = (1.5*tile[k+HALO_WIDTH,l]) / h
		y_out[global_index-1] = (-2.0 * tile[k+HALO_WIDTH,l]) / h
		y_out[global_index-2] = (0.5 * tile[k+HALO_WIDTH,l])/h
    end

    sync_threads()

    nothing

end

# function FACEtoVOL_GPU_shared(y_out, u_face, face, Nx, Ny, ::Val{TILE_DIM}) where {TILE_DIM} #, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
# 	# tidx = threadIdx().x
# 	# i = (blockIdx().x - 1)*TILE_DIM + tidx
# 	i = (blockIdx().x - 1) * TILE_DIM + threadIdx().x
# 	N = Nx*Ny

# 	if i <= N
# 		y_out[i] = 0.0
# 	end

# 	sync_threads()

# 	if face == 1 && i <= N && mod(i,Ny) == 1
# 		y_out[i] = u_face
# 	end

# 	sync_threads()

# 	if face == 2 && i <= N && mod(i,Ny) == 0
# 		y_out[i] = u_face
# 	end

# 	sync_threads()

# 	if face == 3 && i <= Ny
# 		y_out[i] = u_face
# 	end

# 	sync_threads()

# 	if face == 4 && N-Ny+1 <= i <= N
# 		y_out[i] = u_face
# 	end
# 	sync_threads()

# 	nothing
# end

# function VOLtoFACE_GPU_shared(y_in, y_out,  face, Nx, Ny, ::Val{TILE_DIM}) where {TILE_DIM} #, ::Val{TILE_DIM1}, ::Val{TILE_DIM2}) where {TILE_DIM1, TILE_DIM2}
# 	# tidx = threadIdx().x
# 	# i = (blockIdx().x - 1)*TILE_DIM + tidx
# 	i = (blockIdx().x - 1) * TILE_DIM + threadIdx().x
# 	N = Nx*Ny

# 	if i <= N
# 		y_out[i] = 0.0
# 	end

# 	sync_threads()

# 	if face == 1 && i <= N && mod(i,Ny) == 1
# 		y_out[i] = u_face
# 	end

# 	sync_threads()

# 	if face == 2 && i <= N && mod(i,Ny) == 0
# 		y_out[i] = u_face
# 	end

# 	sync_threads()

# 	if face == 3 && i <= Ny
# 		y_out[i] = u_face
# 	end

# 	sync_threads()

# 	if face == 4 && N-Ny+1 <= i <= N
# 		y_out[i] = u_face
# 	end
# 	sync_threads()

# 	nothing
# end

# tester_function : This gives evaluation on CPU, GPU, GPU_shared
function tester_function(f,Nx,TILE_DIM_1,TILE_DIM_2,TILE_DIM)
    Ny = Nx
	@show f
	@eval gpu_function_naive = $(Symbol(f,"_GPU_naive"))
	@eval gpu_function_shared = $(Symbol(f,"_GPU_shared"))
	@show gpu_function_naive
    @show gpu_function_shared
    h = 1/Nx
	# TILE_DIM_1 = 16
	# TILE_DIM_2 = 2

	u = randn(Nx*Ny)
	y_in = CuArray(u)
	y_out = similar(y_in)
	y_out2 = similar(y_in)

	griddim = (div(Nx,TILE_DIM_1) + 1, div(Ny,TILE_DIM_2) + 1)
	blockdim = (TILE_DIM_1,TILE_DIM_2)

	# TILE_DIM = 32
	THREAD_NUM = TILE_DIM
    BLOCK_NUM = div(Nx * Ny,TILE_DIM)+1

	y = f(u,Nx,Ny,h)
	@cuda threads=THREAD_NUM blocks=BLOCK_NUM gpu_function_naive(y_in, y_out, Nx, Ny, h, Val(TILE_DIM))
    @cuda threads=blockdim blocks=griddim gpu_function_shared(y_in, y_out2, Nx, Ny, h, Val(TILE_DIM_1), Val(TILE_DIM_2))
	@show y ≈ Array(y_out)
	@show y ≈ Array(y_out2)
    #@show y - Array(y_out2)

	rep_times = 10

	t_y = time_ns()
	for i in 1:rep_times
		y = f(u,Nx,Ny,h)
	end
	t_y_end = time_ns()
	t1 = t_y_end - t_y

	memsize = length(u) * sizeof(eltype(u))
	@show Float64(t1)
	@printf("CPU Through-put %20.2f\n", 2 * memsize * rep_times / t1)

	println()

	t_y_out = time_ns()
	for i in 1:rep_times
		@cuda threads=THREAD_NUM blocks=BLOCK_NUM gpu_function_naive(y_in, y_out, Nx, Ny, h, Val(TILE_DIM))
		# @cuda threads=THREAD_NUM blocks=BLOCK_NUM D2y_GPU_v2(y_in, y_out, Nx, Ny, h, Val(TILE_DIM))
	end
	synchronize()
	t_y_out_end = time_ns()
	t2 = t_y_out_end - t_y_out
	@show Float64(t2)
	@show Float64(t1)/Float64(t2)
	@printf("GPU Through-put (naive) %20.2f\n", 2 * memsize * rep_times / t2)

	println()

	t_y_out2 = time_ns()
	for i in 1:rep_times
		@cuda threads=blockdim blocks=griddim gpu_function_shared(y_in, y_out2, Nx, Ny, h, Val(TILE_DIM_1), Val(TILE_DIM_2))
	end
	synchronize()
	t_y_out2_end = time_ns()
	t3 = t_y_out2_end - t_y_out2

	@show Float64(t3)
	@show Float64(t1)/Float64(t3)
	@printf("GPU Through-put (shared memory)%20.2f\n", 2 * memsize * rep_times / t3)

end

# tester_function_v2 : This gives evaluation on CPU, GPU_shared with error output
function tester_function_v2(f,Nx,TILE_DIM_1,TILE_DIM_2)
    Ny = Nx
	@show f
	# @eval gpu_function = $(Symbol(f,"_GPU"))
	@eval gpu_function_shared = $(Symbol(f,"_GPU_shared"))
	# @show gpu_function
    @show gpu_function_shared
    h = 1/Nx
	# TILE_DIM_1 = 16
	# TILE_DIM_2 = 2

	u = randn(Nx*Ny)
	y_in = CuArray(u)
	# y_out = similar(y_in)
	y_out2 = similar(y_in)

	griddim = (div(Nx,TILE_DIM_1) + 1, div(Ny,TILE_DIM_2) + 1)
	blockdim = (TILE_DIM_1,TILE_DIM_2)

	# TILE_DIM = 32
	# THREAD_NUM = TILE_DIM
    # BLOCK_NUM = div(Nx * Ny,TILE_DIM)+1

	y = f(u,Nx,Ny,h)
	# @cuda threads=THREAD_NUM blocks=BLOCK_NUM gpu_function(y_in, y_out, Nx, Ny, h, Val(TILE_DIM))
    @cuda threads=blockdim blocks=griddim gpu_function_shared(y_in, y_out2, Nx, Ny, h, Val(TILE_DIM_1), Val(TILE_DIM_2))
	# @show y ≈ Array(y_out)
	@show y ≈ Array(y_out2)

	@show u
	println()
	@show y
	println()
	@show Array(y_out2)
	println()
	@show y - Array(y_out2)

	@show y ≈ Array(y_out2)

	rep_times = 10

	t_y = time_ns()
	for i in 1:rep_times
		y = f(u,Nx,Ny,h)
	end
	t_y_end = time_ns()
	t1 = t_y_end - t_y

	memsize = length(u) * sizeof(eltype(u))
	@show Float64(t1)
	@printf("CPU Through-put %20.2f\n", 2 * memsize * rep_times / t1)

	println()

	# t_y_out = time_ns()
	# for i in 1:rep_times
	# 	@cuda threads=THREAD_NUM blocks=BLOCK_NUM gpu_function(y_in, y_out, Nx, Ny, h, Val(TILE_DIM))
	# 	# @cuda threads=THREAD_NUM blocks=BLOCK_NUM D2y_GPU_v2(y_in, y_out, Nx, Ny, h, Val(TILE_DIM))
	# end
	# synchronize()
	# t_y_out_end = time_ns()
	# t2 = t_y_out_end - t_y_out
	# @show Float64(t2)
	# @show Float64(t1)/Float64(t2)
	# @printf("GPU Through-put (naive) %20.2f\n", 2 * memsize * rep_times / t2)

	# println()

	t_y_out2 = time_ns()
	for i in 1:rep_times
		@cuda threads=blockdim blocks=griddim gpu_function_shared(y_in, y_out2, Nx, Ny, h, Val(TILE_DIM_1), Val(TILE_DIM_2))
	end
	synchronize()
	t_y_out2_end = time_ns()
	t3 = t_y_out2_end - t_y_out2

	@show Float64(t3)
	@show Float64(t1)/Float64(t3)
	@printf("GPU Through-put (shared memory)%20.2f\n", 2 * memsize * rep_times / t3)

end

# tester_function_v3 : This gives evaluation on CPU, GPU_shared with error output turned off
function tester_function_v3(f,Nx,TILE_DIM_1,TILE_DIM_2)
    Ny = Nx
	@show f
	# @eval gpu_function = $(Symbol(f,"_GPU"))
	@eval gpu_function_shared = $(Symbol(f,"_GPU_shared"))
	# @show gpu_function
    @show gpu_function_shared
    h = 1/Nx
	# TILE_DIM_1 = 16
	# TILE_DIM_2 = 2

	u = randn(Nx*Ny)
	y_in = CuArray(u)
	# y_out = similar(y_in)
	y_out2 = similar(y_in)

	griddim = (div(Nx,TILE_DIM_1) + 1, div(Ny,TILE_DIM_2) + 1)
	blockdim = (TILE_DIM_1,TILE_DIM_2)

	@show blockdim
	@show griddim

	# TILE_DIM = 32
	# THREAD_NUM = TILE_DIM
    # BLOCK_NUM = div(Nx * Ny,TILE_DIM)+1

	y = f(u,Nx,Ny,h)
	# @cuda threads=THREAD_NUM blocks=BLOCK_NUM gpu_function(y_in, y_out, Nx, Ny, h, Val(TILE_DIM))
    @cuda threads=blockdim blocks=griddim gpu_function_shared(y_in, y_out2, Nx, Ny, h, Val(TILE_DIM_1), Val(TILE_DIM_2))
	# @show y ≈ Array(y_out)
	@show y ≈ Array(y_out2)

	# @show y - Array(y_out2)

	rep_times = 10

	t_y = time_ns()
	for i in 1:rep_times
		y = f(u,Nx,Ny,h)
	end
	t_y_end = time_ns()
	t1 = t_y_end - t_y

	memsize = length(u) * sizeof(eltype(u))
	@show Float64(t1)
	@printf("CPU Through-put %20.2f\n", 2 * memsize * rep_times / t1)

	println()

	# t_y_out = time_ns()
	# for i in 1:rep_times
	# 	@cuda threads=THREAD_NUM blocks=BLOCK_NUM gpu_function(y_in, y_out, Nx, Ny, h, Val(TILE_DIM))
	# 	# @cuda threads=THREAD_NUM blocks=BLOCK_NUM D2y_GPU_v2(y_in, y_out, Nx, Ny, h, Val(TILE_DIM))
	# end
	# synchronize()
	# t_y_out_end = time_ns()
	# t2 = t_y_out_end - t_y_out
	# @show Float64(t2)
	# @show Float64(t1)/Float64(t2)
	# @printf("GPU Through-put (naive) %20.2f\n", 2 * memsize * rep_times / t2)

	# println()

	t_y_out2 = time_ns()
	for i in 1:rep_times
		@cuda threads=blockdim blocks=griddim gpu_function_shared(y_in, y_out2, Nx, Ny, h, Val(TILE_DIM_1), Val(TILE_DIM_2))
	end
	synchronize()
	t_y_out2_end = time_ns()
	t3 = t_y_out2_end - t_y_out2

	@show Float64(t3)
	@show Float64(t1)/Float64(t3)
	@printf("GPU Through-put (shared memory)%20.2f\n", 2 * memsize * rep_times / t3)

end

function tester_function_FV(f,u_face,face,Nx,TILE_DIM)
	Ny = Nx
	N = Ny*Nx
	@show f
	# @eval gpu_function = $(Symbol(f,"_GPU"))
	@eval gpu_function_shared = $(Symbol(f,"_GPU_shared"))
	# @show gpu_function
    @show gpu_function_shared
	y = zeros(N)
	y_out = CuArray(y)

	# u_face = 2
	# face = 2

	THREAD_NUM = TILE_DIM
	BLOCK_NUM = (N + TILE_DIM-1, TILE_DIM)

	y = f(u_face, face, Nx, Ny)

	@cuda threads=THREAD_NUM blocks=BLOCK_NUM gpu_function_shared(y_out,u_face,face,Nx,Ny,Val(TILE_DIM))

	@show y ≈ Array(y_out)
	# @show y
	# @show Array(y_out)
	# @show y - Array(y_out)

	rep_times = 10
	memsize = length(y) * sizeof(eltype(y))

	t_y = time_ns()
	for i in 1:rep_times
		y .= f(u_face, face, Nx, Ny)
	end
	t_y_end = time_ns()
	t1 = t_y_end - t_y
	@show Float64(t1)
	@printf("CPU Through-put %20.2f\n", 2 * memsize * rep_times / t1)

	t_y_out = time_ns()
	for i in rep_times
		@cuda threads=THREAD_NUM blocks=BLOCK_NUM gpu_function_shared(y_out,u_face,face,Nx,Ny,Val(TILE_DIM))
	end
	synchronize()
	t_y_out_end = time_ns()
	t2 = t_y_out_end - t_y_out
	@show Float64(t2)
	@printf("GPU Through-put (naive) %20.2f\n", 2 * memsize * rep_times / t2)

end

# tester_function(D2x, 1000, 16, 4, 32)
# tester_function(D2y, 1000, 4, 16, 32)

# tester_function_v3(Dx, 10, 16, 4)

# tester_function_v3(Dx, 10, 16, 4)
# tester_function_v3(Dy, 10, 16, 4)

# tester_function_FV(FACEtoVOL,1,1,10,32)

# Nx = 10
# Ny = Nx
# h = 1/Nx
# TILE_DIM_1 = 4
# TILE_DIM_2 = 16

# u = randn(Nx*Ny);
# y_in = CuArray(u);
# y_out = similar(y_in);
# y_out2 = similar(y_in);

# y = D2y(u,Nx,Ny,h)

# griddim = (div(Nx,TILE_DIM_1) + 1, div(Ny,TILE_DIM_2) + 1)
# blockdim = (TILE_DIM_1,TILE_DIM_2)

# # TILE_DIM = 32
# # THREAD_NUM = 32
# # BLOCK_NUM = div(Nx * Ny,TILE_DIM) + 1

# @cuda threads=blockdim blocks=griddim D2y_GPU_v7(y_in, y_out, Nx, Ny, h, Val(TILE_DIM_1), Val(TILE_DIM_2))
# Array(y_out)

# Array(y_out) ≈ y

# D2y(u,Nx,Ny,h)


text:
using Test, Random
Random.seed!(1234321)

using HyperDualMatrixTools
using HyperDualNumbers, LinearAlgebra, SparseArrays, SuiteSparse

@testset "Testing HyperDualMatrixTools" begin
    # Chose a size for matrices
    n = 10

    x = randn(n)                 # real-valued      │
    y = randn(n)                 # real-valued      │
    z = randn(n)                 # real-valued      ├─ vector
    w = randn(n)                 # real-valued      │
    d = x + ε₁*y + ε₂*z + ε₁ε₂*w # hyperdual-valued │

    A = randn(n, n)              # real-valued      │
    B = randn(n, n)              # real-valued      │
    C = randn(n, n)              # real-valued      ├─ matrix
    D = randn(n, n)              # real-valued      │
    M = A + ε₁*B + ε₂*C + ε₁ε₂*D # hyperdual-valued │

    @testset "Testing full matrices" begin
        @testset "Testing `\\` without factorization" begin
            @test A \ (A * x) ≈ x
            @test M \ (M * x) ≈ x
            @test A \ (A * d) ≈ d
            @test M \ (M * d) ≈ d
            @test A * (A \ x) ≈ x
            @test M * (M \ x) ≈ x
            @test A * (A \ d) ≈ d
            @test M * (M \ d) ≈ d
        end

        @testset "Testing `\\` with factorization" begin
            Af = factorize(A)
            Mf = factorize(M)
            @test Af \ (A * x) ≈ x
            @test Mf \ (M * x) ≈ x
            @test Af \ (A * d) ≈ d
            @test Mf \ (M * d) ≈ d
            @test A * (Af \ x) ≈ x
            @test M * (Mf \ x) ≈ x
            @test A * (Af \ d) ≈ d
            @test M * (Mf \ d) ≈ d
        end
    end

    @testset "Inplace factorization" begin
        Mf1 =factorize(M)
        Mf2 = factorize(M)
        factorize(Mf2, 2M)
        @test Mf2.Af == Mf1.Af
        @test Mf2.B == 2Mf1.B
        @test Mf2.C == 2Mf1.C
        @test Mf2.D == 2Mf1.D
        factorize(Mf2, 2M, update_factors=true)
        @test Mf2.Af ≠ Mf1.Af
        @test Mf2.B == 2Mf1.B
        @test Mf2.C == 2Mf1.C
        @test Mf2.D == 2Mf1.D
    end

    @testset "Testing sparse matrices" begin
        spA = sparse(A)  # │
        spB = sparse(B)  # │
        spC = sparse(C)  # ├─ sparse matrices
        spD = sparse(D)  # │
        spM = sparse(M)  # │

        @testset "Check that `ε * A` works" begin
            @test spM ≈ spA + ε₁*spB + ε₂*spC + ε₁ε₂*spD
        end

        @testset "Check that `\\` works without factorization" begin
            @test spA \ (spA * x) ≈ x
            @test spM \ (spM * x) ≈ x
            @test spA \ (spA * d) ≈ d
            @test spM \ (spM * d) ≈ d
            @test spA * (spA \ x) ≈ x
            @test spM * (spM \ x) ≈ x
            @test spA * (spA \ d) ≈ d
            @test spM * (spM \ d) ≈ d
        end

        @testset "Check that `\\` works with factorization" begin
            spAf = factorize(spA)
            spMf = factorize(spM)
            @test spAf \ (spA * x) ≈ x
            @test spMf \ (spM * x) ≈ x
            @test spAf \ (spA * d) ≈ d
            @test spMf \ (spM * d) ≈ d
            @test spA * (spAf \ x) ≈ x
            @test spM * (spMf \ x) ≈ x
            @test spA * (spAf \ d) ≈ d
            @test spM * (spMf \ d) ≈ d
        end
    end

    # TODO Add tests specific to `lu`, `qr`, `cholesky`, etc.

    @testset "Testing $f" for f in [:adjoint, :transpose]
        spA = sparse(A)  # │
        spB = sparse(B)  # │
        spC = sparse(C)  # ├─ sparse matrices
        spD = sparse(D)  # │
        spM = sparse(M)  # │

        # Adjoints / transposes
        x₂ = randn(n)                     # real-valued │
        y₂ = randn(n)                     # real-valued │
        z₂ = randn(n)                     # real-valued ├─ vector
        w₂ = randn(n)                     # real-valued │
        d₂ = x₂ + ε₁*y₂ + ε₂*z₂ + ε₁ε₂*w₂ # dual-valued │

        @eval begin
            x₂′ = $f($x₂)
            x′ = $f($x)
            d₂′ = $f($d₂)
            d′ = $f($d)
            A′  = $f($A)
            M′  = $f($M)
            spA′  = $f($spA)
            spM′  = $f($spM)
        end

        # Check that `\` works with adjoints
        @testset "Check that `\\` works with $f" begin
            @test (x₂′ * (A \ x)) ≈ (x′ * (A′ \ x₂))
            @test (d₂′ * (M \ d)) ≈ (d′ * (M′ \ d₂))
            @test (x₂′ * (spA \ x)) ≈ (x′ * (spA′ \ x₂))
            @test (d₂′ * (spM \ d)) ≈ (d′ * (spM′ \ d₂))
        end

        # Check that `\` works with adjoints of factorized versions
        Af = factorize(A)
        Mf = factorize(M)
        spAf = factorize(spA)
        spMf = factorize(spM)
        @eval begin
            Af′  = $f($Af)
            Mf′  = $f($Mf)
            spAf′  = $f($spAf)
            spMf′  = $f($spMf)
        end
        @testset "Check that factorized version with $f" begin
            @test (x₂′ * (Af \ x)) ≈ (x′ * (Af′ \ x₂))
            @test (d₂′ * (Mf \ d)) ≈ (d′ * (Mf′ \ d₂))
            @test (x₂′ * (spAf \ x)) ≈ (x′ * (spAf′ \ x₂))
            @test (d₂′ * (spMf \ d)) ≈ (d′ * (spMf′ \ d₂))
        end
    end

    @testset "Testing isapprox function" begin
        A2 = A .+ 0ε₁
        @test A2 ≈ A
        @test A ≈ A2
        @test (1.0 + 0ε₁) ≈ 1.0
        @test 1.0 ≈ (1.0 + 0ε₁)

        @test ~(M ≈ M .+ ε₁)
    end
end

println("HyperDualMatrixTools tests finished!")

ext: jl


text:
"""
    struct Cone{T,CO,TR,TP} <: AbstractVolumePrimitive{T, CO}

Volume primitive describing a [Cone](@ref) with its top and bottom circular base
being aligned with the `xy` plane (before possible rotations).

## Parametric types
* `T`: Precision type.
* `CO`: Describes whether the surface belongs to the primitive.
    It can be `ClosedPrimitive`, i.e. the surface points belong to the primitive,
    or `OpenPrimitive`, i.e. the surface points do not belong to the primitive.
* `TR`: Type of the radius `r`.
    * `TR == T`: Cylinder (constant radius `r` at all `z`).
    * `TR == Tuple{T, T}`: Tube (inner radius at `r[1]`, outer radius at `r[2]`).
    * `TR == Tuple{Tuple{T}, Tuple{T}}`: Varying Cylinder (full cylinder with radius changing linearly in `z` from `r[1]` at the bottom to `r[2]` at the top).
    * `TR == Tuple{Tuple{T, T}, Tuple{T, T}}`: Varying Tube (inner radius changes linearly in `z` from `r[1][1]` at the bottom to `r[1][2]` at the top, outer radius changes linearly in `z` from `r[2][1]` at the bottom to `r[1][2]` at the top).
    * `TR == Tuple{Nothing, Tuple{T, T}}`: Cone (Tip at the bottom, top is a circular base with inner radius `r[2][1]` and outer radius `r[2][2]`).
    * `TR == Tuple{Tuple{T, T}, Nothing}`: Cone (Tip at the top, bottom is a circular base with inner radius `r[1][1]` and outer radius `r[1][2]`).
* `TP`: Type of the angular range `φ`.
    * `TP == Nothing`: Full 2π Cone.
    * `TP == T`: Partial Cone ranging from `0` to `φ`.

## Fields
* `r::TR`: Definition of the radius of the `Cone` (in m).
* `φ::TP`: Range in polar angle `φ` over which the `Cone` extends (in radians).
* `hZ::T`: Half of the height of the `Cone` (in m).
* `origin::CartesianPoint{T}`: The position of the center of the `Cone` at the middle height.
* `rotation::SMatrix{3,3,T,9}`: Matrix that describes a rotation of the `Cone` around its `origin`.

## Definition in Configuration File

A `Cone` is defined in the configuration file as part of the `geometry` field
of an object through the field `cone` (or `tube`).

Example definitions of a cylinder looks like this:
```yaml
tube:
  r:
    from: 1.0
    to: 2.0  # => r = (1.0, 2.0)
  h: 2.0     # => hZ = 1.0
  origin:
    z: 1.0   # => origin = [0.0, 0.0, 1.0]
```
This is a hollow cylinder with inner radius 1 at outer radius 2
with a height of 2 and extending over full 2π (no `phi` given),
translated 1 along the `z` axis.

If the radius is not constant over `z`, the `r` entries are divided
into `bottom` and `top`, where `bottom` describes the inner and outer
radius at the bottom circular base and `top` describes the inner and
outer radius at the top circular base (before rotations)
```yaml
cone:
  r:
    bottom:
      from: 1.0
      to: 2.0
    top:
      from: 1.0
      to: 4.0     # => r = ((1.0, 2.0), (1.0, 4.0))
  phi:
    from: 0.0°
    to: 180.0°    # => φ = π
  h: 2.0          # => hZ = 1.0
```
This is half a Cone (`φ` goes from 0 to 180°, i.e. only positive `y` are allowed)
with a height of 2, constant inner radius of 1 and an outer radius which increases
from 2 at the bottom to 4 at the top circular base.

!!! note
    The names `tube` and `cone` in the configuration files are interchangeable.

See also [Constructive Solid Geometry (CSG)](@ref).
"""
struct Cone{T,CO,TR,TP<:Union{Nothing,T}} <: AbstractVolumePrimitive{T, CO}
    r::TR
    φ::TP
    hZ::T

    origin::CartesianPoint{T}
    rotation::SMatrix{3,3,T,9}
end

#Type conversion happens here
function Cone{T,CO}(r, φ, hZ, origin, rotation) where {T,CO}
    _r = _csg_convert_args(T, r)
    (_φ, _rotation) = _handle_phi(_csg_convert_args(T, φ), rotation)
    _hZ = _csg_convert_args(T, hZ)
    Cone{T,CO,typeof(_r),typeof(_φ)}(_r, _φ, _hZ, origin, _rotation)
end

#Type promotion happens here
function Cone(CO, r::TR, φ::TP, hZ::TZ, origin::PT, rotation::ROT) where {TR, TP, TZ, PT, ROT}
    eltypes = _csg_get_promoted_eltype.((TR, TZ, TP, PT, ROT))
    T = float(promote_type(eltypes...))
    Cone{T,CO}(r, φ, hZ, origin, rotation)
end

function Cone(::Type{CO}=ClosedPrimitive;
    # define default parameters as Int to not influence type promotion
    r = 1,
    φ = nothing,
    hZ = 1,
    origin = zero(CartesianPoint{Int}),
    rotation = one(SMatrix{3, 3, Int, 9})
) where {CO}
    Cone(CO, r, φ, hZ, origin, rotation)
end

function Cone{T}(::Type{CO}=ClosedPrimitive;
    r = 1.0,
    φ = nothing,
    hZ = 1.0,
    origin = zero(CartesianPoint{Float64}),
    rotation = one(SMatrix{3, 3, Float64, 9})
) where {T, CO}
    Cone{T,CO}(r, φ, hZ, origin, rotation)
end

Cone{T,CO,TR,TP}( c::Cone{T,CO,TR,TP}; COT = CO,
            origin::CartesianPoint{T} = c.origin,
            rotation::SMatrix{3,3,T,9} = c.rotation) where {T,CO<:Union{ClosedPrimitive, OpenPrimitive},TR,TP<:Union{Nothing,T}} =
    Cone{T,COT,TR,TP}(c.r, c.φ, c.hZ, origin, rotation)

####################################################################
####################################################################

#=
___
   |
   |
   |
___|

=#
const Cylinder{T,CO} = Cone{T,CO,T,Nothing} # Full in φ
const PartialCylinder{T,CO} = Cone{T,CO,T,T}

### Cylinder

function _in(pt::CartesianPoint, c::Cylinder{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    az = abs(pt.z)
    az <= c.hZ + csgtol && begin
        r = hypot(pt.x, pt.y)
        r <= c.r + csgtol
    end
end
function _in(pt::CartesianPoint, c::Cylinder{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.z) < c.hZ - csgtol &&
    hypot(pt.x, pt.y) < c.r - csgtol
end

function surfaces(t::Cylinder{T}) where {T}
    bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
    top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
    mantle = FullConeMantle{T,:inwards}((t.r, t.r), t.φ, t.hZ, t.origin, t.rotation)
    e_bot = CircularArea{T}(r = t.r, φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
    e_top = CircularArea{T}(r = t.r, φ = t.φ, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
    # normals of the surfaces show inside the volume primitives.
    e_top, e_bot, mantle
end

### PartialCylinder

function _in(pt::CartesianPoint, c::PartialCylinder{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    az = abs(pt.z)
    az <= c.hZ + csgtol && begin
        r = hypot(pt.x, pt.y)
        r <= c.r + csgtol &&
        _in_angular_interval_closed(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)
    end
end
function _in(pt::CartesianPoint, c::PartialCylinder{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.z) + csgtol < c.hZ && begin
        r = hypot(pt.x, pt.y)
        csgtol + r < c.r &&
        _in_angular_interval_open(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)
    end
end

function surfaces(t::PartialCylinder{T}) where {T}
    bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
    top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
    mantle = PartialConeMantle{T,:inwards}((t.r, t.r), t.φ, t.hZ, t.origin, t.rotation)
    e_bot = PartialCircularArea{T}(r = t.r, φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
    e_top = PartialCircularArea{T}(r = t.r, φ = t.φ, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
    poly_l = _transform_into_global_coordinate_system(Quadrangle{T}((
        CartesianPoint(CylindricalPoint{T}(zero(T), zero(T), -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(zero(T), zero(T), +t.hZ)),
        CartesianPoint(CylindricalPoint{T}(    t.r, zero(T), +t.hZ)),
        CartesianPoint(CylindricalPoint{T}(    t.r, zero(T), -t.hZ)) )), t)
    poly_r = _transform_into_global_coordinate_system(Quadrangle{T}((
        CartesianPoint(CylindricalPoint{T}(zero(T), t.φ, -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(    t.r, t.φ, -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(    t.r, t.φ, +t.hZ)),
        CartesianPoint(CylindricalPoint{T}(zero(T), t.φ, +t.hZ)) )), t)
    # normals of the surfaces show inside the volume primitives.
    (e_top, e_bot, mantle, poly_l, poly_r)
end

####################################################################
####################################################################

#=
___
   \
    \
_____\

=#
const VaryingCylinder{T,CO} = Cone{T,CO,Tuple{Tuple{T},Tuple{T}},Nothing} # Full in φ
const PartialVaryingCylinder{T,CO} = Cone{T,CO,Tuple{Tuple{T},Tuple{T}},T}

function _in(pt::CartesianPoint, c::VaryingCylinder{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    az = abs(pt.z)
    az <= c.hZ + csgtol && begin
        r = hypot(pt.x, pt.y)
        rz = radius_at_z(c.hZ, c.r[1][1], c.r[2][1], pt.z)
        r <= rz + csgtol
    end
end
function _in(pt::CartesianPoint, c::VaryingCylinder{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.z) + csgtol < c.hZ &&
    csgtol + hypot(pt.x, pt.y) < radius_at_z(c.hZ, c.r[1][1], c.r[2][1], pt.z)
end

function surfaces(t::VaryingCylinder{T}) where {T}
    bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
    top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
    mantle = FullConeMantle{T,:inwards}((t.r[1][1], t.r[2][1]), t.φ, t.hZ, t.origin, t.rotation)
    e_bot = CircularArea{T}(r = t.r[1][1], φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
    e_top = CircularArea{T}(r = t.r[2][1], φ = t.φ, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
    # normals of the surfaces show inside the volume primitives.
    e_top, e_bot, mantle
end

function _in(pt::CartesianPoint, c::PartialVaryingCylinder{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    az = abs(pt.z)
    az <= c.hZ + csgtol && begin
        r = hypot(pt.x, pt.y)
        rz = radius_at_z(c.hZ, c.r[1][1], c.r[2][1], pt.z)
        r <= rz + csgtol &&
        _in_angular_interval_closed(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)
    end
end
function _in(pt::CartesianPoint, c::PartialVaryingCylinder{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.z) + csgtol < c.hZ && begin
        r = hypot(pt.x, pt.y)
        csgtol + r < radius_at_z(c.hZ, c.r[1][1], c.r[2][1], pt.z) &&
        _in_angular_interval_open(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)
    end
end

function surfaces(t::PartialVaryingCylinder{T}) where {T}
    bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
    top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
    mantle = PartialConeMantle{T,:inwards}((t.r[1][1], t.r[2][1]), t.φ, t.hZ, t.origin, t.rotation)
    e_bot = PartialCircularArea{T}(r = t.r[1][1], φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
    e_top = PartialCircularArea{T}(r = t.r[2][1], φ = t.φ, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
    poly_l = _transform_into_global_coordinate_system(Quadrangle{T}((
        CartesianPoint(CylindricalPoint{T}(zero(T),    zero(T), -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(zero(T),    zero(T), +t.hZ)),
        CartesianPoint(CylindricalPoint{T}( t.r[2][1], zero(T), +t.hZ)),
        CartesianPoint(CylindricalPoint{T}( t.r[1][1], zero(T), -t.hZ)) )), t)
    poly_r = _transform_into_global_coordinate_system(Quadrangle{T}((
        CartesianPoint(CylindricalPoint{T}(zero(T),    t.φ, -t.hZ)),
        CartesianPoint(CylindricalPoint{T}( t.r[1][1], t.φ, -t.hZ)),
        CartesianPoint(CylindricalPoint{T}( t.r[2][1], t.φ, +t.hZ)),
        CartesianPoint(CylindricalPoint{T}(zero(T),    t.φ, +t.hZ)) )), t)
    # normals of the surfaces show inside the volume primitives.
    e_top, e_bot, mantle, poly_l, poly_r
end

####################################################################
####################################################################

#=
  ___
 /   \
/_____\

=#
const VaryingTube{T,CO} = Cone{T,CO,Tuple{Tuple{T,T},Tuple{T,T}},Nothing} # Full in φ
const PartialVaryingTube{T,CO} = Cone{T,CO,Tuple{Tuple{T,T},Tuple{T,T}},T}

#(r_bot_in = r[1][1], r_bot_out = r[1][2], r_top_in = r[2][1], r_top_out = r[2][2])

function _in(pt::CartesianPoint, c::VaryingTube{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    az = abs(pt.z)
    az <= c.hZ + csgtol && begin
        r = hypot(pt.x, pt.y)
        r_in = radius_at_z(c.hZ, c.r[1][1], c.r[2][1], pt.z)
        r_out = radius_at_z(c.hZ, c.r[1][2], c.r[2][2], pt.z)
        r_in - csgtol <= r  &&
        r <= r_out + csgtol
    end
end
function _in(pt::CartesianPoint, c::VaryingTube{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.z) + csgtol < c.hZ &&
    csgtol + radius_at_z(c.hZ, c.r[1][1], c.r[2][1], pt.z) < hypot(pt.x, pt.y) < radius_at_z(c.hZ, c.r[1][2], c.r[2][2], pt.z) - csgtol
end

function surfaces(t::VaryingTube{T}) where {T}
    bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
    top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
    inner_mantle = FullConeMantle{T,:outwards}((t.r[1][1], t.r[2][1]), t.φ, t.hZ, t.origin, t.rotation)
    outer_mantle = FullConeMantle{T,:inwards}( (t.r[1][2], t.r[2][2]), t.φ, t.hZ, t.origin, t.rotation)
    e_bot = Annulus{T}(r = (t.r[1][1], t.r[1][2]), φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
    e_top = Annulus{T}(r = (t.r[2][1], t.r[2][2]), φ = t.φ, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
    # normals of the surfaces show inside the volume primitives.
    e_top, e_bot, inner_mantle, outer_mantle
end

function _in(pt::CartesianPoint, c::PartialVaryingTube{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    az = abs(pt.z)
    az <= c.hZ + csgtol && begin
        r = hypot(pt.x, pt.y)
        r_in = radius_at_z(c.hZ, c.r[1][1], c.r[2][1], pt.z)
        r_out = radius_at_z(c.hZ, c.r[1][2], c.r[2][2], pt.z)
        r_in - csgtol <= r &&
        r <= r_out + csgtol &&
        _in_angular_interval_closed(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)
    end
end
function _in(pt::CartesianPoint, c::PartialVaryingTube{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.z) + csgtol < c.hZ && begin
        r = hypot(pt.x, pt.y)
        csgtol + radius_at_z(c.hZ, c.r[1][1], c.r[2][1], pt.z) < r < radius_at_z(c.hZ, c.r[1][2], c.r[2][2], pt.z) - csgtol &&
        _in_angular_interval_open(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)
    end
end

function surfaces(t::PartialVaryingTube{T}) where {T}
    bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
    top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
    inner_mantle = PartialConeMantle{T,:outwards}((t.r[1][1], t.r[2][1]), t.φ, t.hZ, t.origin, t.rotation)
    outer_mantle = PartialConeMantle{T,:inwards}( (t.r[1][2], t.r[2][2]), t.φ, t.hZ, t.origin, t.rotation)
    e_bot = PartialAnnulus{T}(r = (t.r[1][1], t.r[1][2]), φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
    e_top = PartialAnnulus{T}(r = (t.r[2][1], t.r[2][2]), φ = t.φ, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
    poly_l = _transform_into_global_coordinate_system(Quadrangle{T}((
        CartesianPoint(CylindricalPoint{T}(t.r[1][1], zero(T), -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(t.r[2][1], zero(T), +t.hZ)),
        CartesianPoint(CylindricalPoint{T}(t.r[2][2], zero(T), +t.hZ)),
        CartesianPoint(CylindricalPoint{T}(t.r[1][2], zero(T), -t.hZ)) )), t)
    poly_r = _transform_into_global_coordinate_system(Quadrangle{T}((
        CartesianPoint(CylindricalPoint{T}(t.r[1][1], t.φ, -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(t.r[1][2], t.φ, -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(t.r[2][2], t.φ, +t.hZ)),
        CartesianPoint(CylindricalPoint{T}(t.r[2][1], t.φ, +t.hZ)) )), t)
    # normals of the surfaces show inside the volume primitives.
    e_top, e_bot, inner_mantle, outer_mantle, poly_l, poly_r
end

####################################################################
####################################################################

#=

\
 \
__\

=#
const UpwardCone{T,CO} = Cone{T,CO,Tuple{Tuple{Nothing,T},Nothing},Nothing} # Full in φ
const PartialUpwardCone{T,CO} = Cone{T,CO,Tuple{Tuple{Nothing,T},Nothing},T}

function _in(pt::CartesianPoint, c::UpwardCone{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    az = abs(pt.z)
    az <= c.hZ + csgtol && begin
        r = hypot(pt.x, pt.y)
        r_out = radius_at_z(c.hZ, c.r[1][2], zero(T), pt.z)
        r <= r_out + csgtol
    end
end
function _in(pt::CartesianPoint, c::UpwardCone{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.z) + csgtol < c.hZ &&
    hypot(pt.x, pt.y) + csgtol < radius_at_z(c.hZ, c.r[1][2], zero(T), pt.z)
end

function surfaces(t::UpwardCone{T}) where {T}
    bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
    outer_mantle = FullConeMantle{T,:inwards}( (t.r[1][2], zero(T)), t.φ, t.hZ, t.origin, t.rotation)
    e_bot = CircularArea{T}(r = t.r[1][2], φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
    # normals of the surfaces show inside the volume primitives.
    e_bot, outer_mantle
end

function _in(pt::CartesianPoint, c::PartialUpwardCone{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    az = abs(pt.z)
    az <= c.hZ + csgtol && begin
        r = hypot(pt.x, pt.y)
        r_out = radius_at_z(c.hZ, c.r[1][2], zero(T), pt.z)
        r <= r_out + csgtol &&
        _in_angular_interval_closed(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)
    end
end
function _in(pt::CartesianPoint, c::PartialUpwardCone{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.z) + csgtol < c.hZ && begin
        r = hypot(pt.x, pt.y)
        csgtol + r < radius_at_z(c.hZ, c.r[1][2], zero(T), pt.z) &&
        _in_angular_interval_open(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)
    end
end

function surfaces(t::PartialUpwardCone{T}) where {T}
    bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
    outer_mantle = PartialConeMantle{T,:inwards}( (t.r[1][2], zero(T)), t.φ, t.hZ, t.origin, t.rotation)
    e_bot = PartialCircularArea{T}(r = t.r[1][2], φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
    poly_l = _transform_into_global_coordinate_system(Triangle{T}((
        CartesianPoint(CylindricalPoint{T}(zero(T),   zero(T), -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(zero(T),   zero(T), +t.hZ)),
        CartesianPoint(CylindricalPoint{T}(t.r[1][2], zero(T), -t.hZ)) )), t)
    poly_r = _transform_into_global_coordinate_system(Triangle{T}((
        CartesianPoint(CylindricalPoint{T}(zero(T),   t.φ, -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(t.r[1][2], t.φ, -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(zero(T),   t.φ, +t.hZ)) )), t)
    # normals of the surfaces show inside the volume primitives.
    e_bot, outer_mantle, poly_l, poly_r
end

####################################################################
####################################################################

#=
___
  /
 /
/

=#
const DownwardCone{T,CO} = Cone{T,CO,Tuple{Nothing,Tuple{Nothing,T}},Nothing} # Full in φ
const PartialDownwardCone{T,CO} = Cone{T,CO,Tuple{Nothing,Tuple{Nothing,T}},T}

function _in(pt::CartesianPoint, c::DownwardCone{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    az = abs(pt.z)
    az <= c.hZ + csgtol && begin
        r = hypot(pt.x, pt.y)
        r_out = radius_at_z(c.hZ, zero(T), c.r[2][2], pt.z)
        r <= r_out + csgtol
    end
end
function _in(pt::CartesianPoint, c::DownwardCone{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.z) + csgtol < c.hZ &&
    hypot(pt.x, pt.y) + csgtol < radius_at_z(c.hZ, zero(T), c.r[2][2], pt.z)
end

function surfaces(t::DownwardCone{T}) where {T}
    top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
    outer_mantle = FullConeMantle{T,:inwards}( (zero(T), t.r[2][2]), t.φ, t.hZ, t.origin, t.rotation)
    e_top = CircularArea{T}(r = t.r[2][2], φ = t.φ, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
    # normals of the surfaces show inside the volume primitives.
    e_top, outer_mantle
end

function _in(pt::CartesianPoint, c::PartialDownwardCone{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    az = abs(pt.z)
    az <= c.hZ + csgtol && begin
        r = hypot(pt.x, pt.y)
        r_out = radius_at_z(c.hZ, zero(T), c.r[2][2], pt.z)
        r <= r_out + csgtol &&
        _in_angular_interval_closed(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)
    end
end
function _in(pt::CartesianPoint, c::PartialDownwardCone{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.z) + csgtol < c.hZ && begin
        r = hypot(pt.x, pt.y)
        csgtol + r < radius_at_z(c.hZ, zero(T), c.r[2][2], pt.z) &&
        _in_angular_interval_open(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)
    end
end

function surfaces(t::PartialDownwardCone{T}) where {T}
    bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
    outer_mantle = PartialConeMantle{T,:inwards}( (t.r[2][2], zero(T)), t.φ, t.hZ, t.origin, t.rotation)
    e_bot = PartialCircularArea{T}(r = t.r[2][2], φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
    poly_l = _transform_into_global_coordinate_system(Triangle{T}((
        CartesianPoint(CylindricalPoint{T}(zero(T),   zero(T), -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(zero(T),   zero(T), +t.hZ)),
        CartesianPoint(CylindricalPoint{T}(t.r[2][2], zero(T), -t.hZ)) )), t)
    poly_r = _transform_into_global_coordinate_system(Triangle{T}((
        CartesianPoint(CylindricalPoint{T}(zero(T),   t.φ, -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(t.r[2][2], t.φ, -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(zero(T),   t.φ, +t.hZ)) )), t)
    # normals of the surfaces show inside the volume primitives.
    e_bot, outer_mantle, poly_l, poly_r
end

####################################################################
####################################################################

#=

  /\
 /  \
/____\

=#
const TopClosedTube{T,CO} = Cone{T,CO,Tuple{Tuple{T,T},T},Nothing} # Full in φ
const PartialTopClosedTube{T,CO} = Cone{T,CO,Tuple{Tuple{T,T},T},T}

function _in(pt::CartesianPoint, c::TopClosedTube{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    az = abs(pt.z)
    az <= c.hZ + csgtol && begin
        r = hypot(pt.x, pt.y)
        r_in = radius_at_z(c.hZ, c.r[1][1], c.r[2], pt.z)
        r_out = radius_at_z(c.hZ, c.r[1][2], c.r[2], pt.z)
        r_in - csgtol <= r &&
        r <= r_out + csgtol
    end
end
function _in(pt::CartesianPoint, c::TopClosedTube{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.z) + csgtol < c.hZ &&
    csgtol + radius_at_z(c.hZ, c.r[1][1], c.r[2], pt.z) < hypot(pt.x, pt.y) < radius_at_z(c.hZ, c.r[1][2], c.r[2], pt.z) - csgtol
end

function surfaces(t::TopClosedTube{T}) where {T}
    bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
    # top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
    inner_mantle = FullConeMantle{T,:outwards}((t.r[1][1], t.r[2]), t.φ, t.hZ, t.origin, t.rotation)
    outer_mantle = FullConeMantle{T,:inwards}( (t.r[1][2], t.r[2]), t.φ, t.hZ, t.origin, t.rotation)
    e_bot = Annulus{T}(r = (t.r[1][1], t.r[1][2]), φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
    # e_top = Annulus{T}(r = (t.r[2][1], t.r[2][2]), φ = t.φ, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
    # normals of the surfaces show inside the volume primitives.
    e_bot, inner_mantle, outer_mantle
end

function _in(pt::CartesianPoint, c::PartialTopClosedTube{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    az = abs(pt.z)
    az <= c.hZ + csgtol && begin
        r = hypot(pt.x, pt.y)
        r_in = radius_at_z(c.hZ, c.r[1][1], c.r[2], pt.z)
        r_out = radius_at_z(c.hZ, c.r[1][2], c.r[2], pt.z)
        r_in - csgtol <= r &&
        r <= r_out + csgtol &&
        _in_angular_interval_closed(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)
    end
end
function _in(pt::CartesianPoint, c::PartialTopClosedTube{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.z) + csgtol < c.hZ && begin
        r = hypot(pt.x, pt.y)
        csgtol + radius_at_z(c.hZ, c.r[1][1], c.r[2], pt.z) < r < radius_at_z(c.hZ, c.r[1][2], c.r[2], pt.z) - csgtol &&
        _in_angular_interval_open(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)
    end
end

function surfaces(t::PartialTopClosedTube{T}) where {T}
    bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
    # top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
    inner_mantle = PartialConeMantle{T,:outwards}((t.r[1][1], t.r[2]), t.φ, t.hZ, t.origin, t.rotation)
    outer_mantle = PartialConeMantle{T,:inwards}( (t.r[1][2], t.r[2]), t.φ, t.hZ, t.origin, t.rotation)
    e_bot = PartialAnnulus{T}(r = (t.r[1][1], t.r[1][2]), φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
    # e_top = Annulus{T}(r = (t.r[2][1], t.r[2][2]), φ = t.φ, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
    poly_l = _transform_into_global_coordinate_system(Triangle{T}((
        CartesianPoint(CylindricalPoint{T}(t.r[1][1], zero(T), -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(   t.r[2], zero(T), +t.hZ)),
        CartesianPoint(CylindricalPoint{T}(t.r[1][2], zero(T), -t.hZ)) )), t)
    poly_r = _transform_into_global_coordinate_system(Triangle{T}((
        CartesianPoint(CylindricalPoint{T}(t.r[1][1], t.φ, -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(t.r[1][2], t.φ, -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(   t.r[2], t.φ, +t.hZ)) )), t)
    # normals of the surfaces show inside the volume primitives.
    e_bot, inner_mantle, outer_mantle, poly_l, poly_r
end

####################################################################
####################################################################

#=

______
\    /
 \  /
  \/

=#
const BottomClosedTube{T,CO} = Cone{T,CO,Tuple{T,Tuple{T,T}},Nothing} # Full in φ
const PartialBottomClosedTube{T,CO} = Cone{T,CO,Tuple{T,Tuple{T,T}},T}

function _in(pt::CartesianPoint, c::BottomClosedTube{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    az = abs(pt.z)
    az <= c.hZ + csgtol && begin
        r = hypot(pt.x, pt.y)
        r_in = radius_at_z(c.hZ, c.r[1], c.r[2][1], pt.z)
        r_out = radius_at_z(c.hZ, c.r[1], c.r[2][2], pt.z)
        r_in - csgtol <= r &&
        r <= r_out + csgtol
    end
end
function _in(pt::CartesianPoint, c::BottomClosedTube{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.z) + csgtol < c.hZ &&
    csgtol + radius_at_z(c.hZ, c.r[1], c.r[2][1], pt.z) < hypot(pt.x, pt.y) < radius_at_z(c.hZ, c.r[1], c.r[2][2], pt.z) - csgtol
end

function surfaces(t::BottomClosedTube{T}) where {T}
    # bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
    top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
    inner_mantle = FullConeMantle{T,:outwards}((t.r[1], t.r[2][1]), t.φ, t.hZ, t.origin, t.rotation)
    outer_mantle = FullConeMantle{T,:inwards}( (t.r[1], t.r[2][2]), t.φ, t.hZ, t.origin, t.rotation)
    # e_bot = Annulus{T}(r = (t.r[1][1], t.r[1][2]), φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
    e_top = Annulus{T}(r = (t.r[2][1], t.r[2][2]), φ = t.φ, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
    # normals of the surfaces show inside the volume primitives.
    e_top, inner_mantle, outer_mantle
end

function _in(pt::CartesianPoint, c::PartialBottomClosedTube{T,ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    az = abs(pt.z)
    az <= c.hZ + csgtol && begin
        r = hypot(pt.x, pt.y)
        r_in = radius_at_z(c.hZ, c.r[1], c.r[2][1], pt.z)
        r_out = radius_at_z(c.hZ, c.r[1], c.r[2][2], pt.z)
        r_in - csgtol <= r &&
        r <= r_out + csgtol &&
        _in_angular_interval_closed(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)
    end
end
function _in(pt::CartesianPoint, c::PartialBottomClosedTube{T,OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T}
    abs(pt.z) + csgtol < c.hZ && begin
        r = hypot(pt.x, pt.y)
        csgtol + radius_at_z(c.hZ, c.r[1], c.r[2][1], pt.z) < r < radius_at_z(c.hZ, c.r[1], c.r[2][2], pt.z) - csgtol &&
        _in_angular_interval_open(atan(pt.y, pt.x), c.φ, csgtol = csgtol / r)
    end
end

function surfaces(t::PartialBottomClosedTube{T}) where {T}
    # bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
    top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
    inner_mantle = PartialConeMantle{T,:outwards}((t.r[1], t.r[2][1]), t.φ, t.hZ, t.origin, t.rotation)
    outer_mantle = PartialConeMantle{T,:inwards}( (t.r[1], t.r[2][2]), t.φ, t.hZ, t.origin, t.rotation)
    # e_bot = PartialAnnulus{T}(r = (t.r[1][1], t.r[1][2]), φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
    e_top = PartialAnnulus{T}(r = (t.r[2][1], t.r[2][2]), φ = t.φ, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
    poly_l = _transform_into_global_coordinate_system(Triangle{T}((
        CartesianPoint(CylindricalPoint{T}(t.r[2][1], zero(T), +t.hZ)),
        CartesianPoint(CylindricalPoint{T}(t.r[2][2], zero(T), +t.hZ)),
        CartesianPoint(CylindricalPoint{T}(   t.r[1], zero(T), -t.hZ)) )), t)
    poly_r = _transform_into_global_coordinate_system(Triangle{T}((
        CartesianPoint(CylindricalPoint{T}(t.r[2][1], t.φ, +t.hZ)),
        CartesianPoint(CylindricalPoint{T}(   t.r[1], t.φ, -t.hZ)),
        CartesianPoint(CylindricalPoint{T}(t.r[2][2], t.φ, +t.hZ)) )), t)
    # normals of the surfaces show inside the volume primitives.
    e_top, inner_mantle, outer_mantle, poly_l, poly_r
end

####################################################################
####################################################################

function Geometry(::Type{T}, t::Type{Cone}, dict::AbstractDict, input_units::NamedTuple, transformations::Transformations{T}) where {T}
    length_unit = input_units.length
    angle_unit = input_units.angle
    origin = get_origin(T, dict, length_unit)
    rotation = get_rotation(T, dict, angle_unit)

    r = parse_r_of_primitive(T, dict, length_unit, Cone)
    φ_interval = parse_φ_of_primitive(T, dict, angle_unit)
    φ = if φ_interval isa Tuple{T,T}
        rotation = rotation * RotZ{T}(φ_interval[1])
        φ_interval[2] - φ_interval[1]
    elseif isnothing(φ_interval)
        nothing
    else
        throw(ConfigFileError("Error when trying to parse φ from configuration file."))
    end

    hZ = if haskey(dict, "h")
        _parse_value(T, dict["h"], length_unit) / 2
    elseif haskey(dict, "z")
        z = parse_height_of_primitive(T, dict, length_unit)
        hZ = typeof(z) <: Real ? z : (z[2] - z[1])/2
        if z isa Real
            @warn "Deprecation warning: Field `z` for `Cone` is deprecated.
                Use `h` instead to specify the height of the primitive."
        else # z isa Tuple
            @warn "Deprecation warning: Field `z` for `Cone` is deprecated.
                Use `h` instead to specify the height of the primitive.
                There might be a conflict with the possible field `origin`:
                The `z` component of the origin of the primitive is overwritten by the `z`."
            origin = CartesianPoint{T}(origin[1], origin[2], mean(z))
       end
       hZ
    end

    cone = Cone{T}(ClosedPrimitive;
        r = r,
        φ = φ,
        hZ = hZ,
        origin = origin,
        rotation = rotation
    )
    return transform(cone, transformations)
end

function Dictionary(c::Cone{T, <:Any, TR})::OrderedDict{String, Any} where {T, TR}
    dict = OrderedDict{String, Any}()
    name  = "tube"
    if TR <: Real
        dict["r"] = c.r
    elseif TR <: Tuple{<:Real, <:Real}
        dict["r"] = OrderedDict{String, Any}("from" => c.r[1], "to" => c.r[2])
    else
        name = "cone"
        dict["r"] = OrderedDict{String, Any}()
        dict["r"]["bottom"] = if c.r[1] isa Tuple{<:Real}
            c.r[1][1]
        elseif isnothing(c.r[1])
            0
        else
            OrderedDict{String, Any}("from" => c.r[1][1], "to" => c.r[1][2])
        end
        dict["r"]["top"] = if c.r[2] isa Tuple{<:Real}
            c.r[2][1]
        elseif isnothing(c.r[2])
            0
        else
            OrderedDict{String, Any}("from" => c.r[2][1], "to" => c.r[2][2])
        end
    end
    if !isnothing(c.φ) dict["phi"] = OrderedDict("from" => "0°", "to" => string(rad2deg(c.φ))*"°") end
    dict["h"] = 2*c.hZ
    if c.origin != zero(CartesianVector{T}) dict["origin"] = c.origin end
    if c.rotation != one(SMatrix{3,3,T,9})
        d = Dictionary(c.rotation)
        if unique(keys(d)) == ["Z"]
            φ0 = mod2pi(_parse_value(T, d["Z"], internal_angle_unit))
            dict["phi"] = OrderedDict("from" => string(rad2deg(φ0))*"°", "to" => string(rad2deg(φ0 + c.φ))*"°")
        else
            dict["rotation"] = d
        end
    end
    OrderedDict{String, Any}(name => dict)
end

extremum(c::Cone{T}) where {T} = hypot(c.hZ, max((c.r...)...))

# Cylinder

# PartialCylinder

# Tube
# function _in(pt::CartesianPoint, c::Tube{T,ClosedPrimitive}) where {T}
#     abs(pt.z) <= c.hZ &&
#     c.r[1] <= hypot(pt.x, pt.y) <= c.r[2]
# end
# function _in(pt::CartesianPoint, c::Tube{T,OpenPrimitive}) where {T}
#     abs(pt.z) < c.hZ &&
#     c.r[1] < hypot(pt.x, pt.y) < c.r[2]
# end

# function surfaces(t::Tube{T}) where {T}
#     bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
#     top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
#     inner_mantle = ConeMantle{T}((t.r[1], t.r[1]), t.φ, t.hZ, t.origin, t.rotation)
#     outer_mantle = ConeMantle{T}((t.r[2], t.r[2]), t.φ, t.hZ, t.origin, t.rotation)
#     e_bot = Annulus{T}(r = t.r, φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
#     e_top = Annulus{T}(r = t.r, φ = t.φ, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
#     # normals of the surfaces show inside the volume primitives.
#     e_top, e_bot, inner_mantle, outer_mantle
# end

# PartialTube
# function _in(pt::CartesianPoint, c::PartialTube{T,ClosedPrimitive}) where {T}
#     abs(pt.z) <= c.hZ &&
#     c.r[1] <= hypot(pt.x, pt.y) <= c.r[2] &&
#     _in_angular_interval_closed(atan(pt.y, pt.x), c.φ)
# end
# function _in(pt::CartesianPoint, c::PartialTube{T,OpenPrimitive}) where {T}
#     abs(pt.z) < c.hZ &&
#     c.r[1] < hypot(pt.x, pt.y) < c.r &&
#     _in_angular_interval_open(atan(pt.y, pt.x), c.φ)
# end

# function surfaces(t::PartialTube{T}) where {T}
#     bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
#     top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
#     inner_mantle = PartialConeMantle{T}((t.r[1], t.r[1]), t.φ, t.hZ, t.origin, t.rotation)
#     outer_mantle = PartialConeMantle{T}((t.r[2], t.r[2]), t.φ, t.hZ, t.origin, t.rotation)
#     e_bot = PartialAnnulus{T}(r = t.r, φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
#     e_top = PartialAnnulus{T}(r = t.r, φ = t.φ, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
#     poly_l = _transform_into_global_coordinate_system(Quadrangle{T}((
#         CartesianPoint(CylindricalPoint{T}( t.r[1], t.φ[1], -t.hZ)),
#         CartesianPoint(CylindricalPoint{T}( t.r[2], t.φ[1], -t.hZ)),
#         CartesianPoint(CylindricalPoint{T}( t.r[2], t.φ[1], +t.hZ)),
#         CartesianPoint(CylindricalPoint{T}( t.r[1], t.φ[1], +t.hZ)) )), t)
#     poly_r = _transform_into_global_coordinate_system(Quadrangle{T}((
#         CartesianPoint(CylindricalPoint{T}( t.r[2], t.φ[2], -t.hZ)),
#         CartesianPoint(CylindricalPoint{T}( t.r[2], t.φ[2], +t.hZ)),
#         CartesianPoint(CylindricalPoint{T}( t.r[1], t.φ[2], +t.hZ)),
#         CartesianPoint(CylindricalPoint{T}( t.r[1], t.φ[2], -t.hZ)) )), t)
#     # normals of the surfaces show inside the volume primitives.
#     e_top, e_bot, inner_mantle, outer_mantle, poly_l, poly_r
# end

# VaryingTube
# function _in(pt::CartesianPoint, c::VaryingTube{T,ClosedPrimitive}) where {T}
#     abs(pt.z) <= c.hZ &&
#     c.r[1] <= hypot(pt.x, pt.y) <= c.r[2]
# end
# function _in(pt::CartesianPoint, c::VaryingTube{T,OpenPrimitive}) where {T}
#     abs(pt.z) < c.hZ &&
#     c.r[1] < hypot(pt.x, pt.y) < c.r[2]
# end

# function surfaces(t::VaryingTube{T}) where {T}
#     bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
#     top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
#     inner_mantle = CylinderMantle{T}(t.r[1], t.φ, t.hZ, t.origin, t.rotation)
#     outer_mantle = CylinderMantle{T}(t.r[2], t.φ, t.hZ, t.origin, t.rotation)
#     e_bot = Annulus{T}(r = t.r, φ = t.φ, origin = bot_center_pt, rotation = t.rotation)
#     e_top = Annulus{T}(r = t.r, φ = t.φ, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
#     # normals of the surfaces show inside the volume primitives.
#     e_top, e_bot, inner_mantle, outer_mantle
# end

# function surfaces(t::VaryingTube{T,CO}) where {T,CO}
#     bot_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), -t.hZ), t)
#     top_center_pt = _transform_into_global_coordinate_system(CartesianPoint{T}(zero(T), zero(T), +t.hZ), t)
#     in_mantle  = ConeMantle{T,Tuple{T,T},Nothing}((t.r[1][1], t.r[2][1]), t.φ, t.hZ, t.origin, t.rotation)
#     out_mantle = ConeMantle{T,Tuple{T,T},Nothing}((t.r[1][2], t.r[2][2]), t.φ, t.hZ, t.origin, t.rotation)
#     e_bot = EllipticalSurface{T,Tuple{T,T},Nothing}(r = t.r[1], φ = nothing, origin = bot_center_pt, rotation = t.rotation)
#     e_top = EllipticalSurface{T,Tuple{T,T},Nothing}(r = t.r[2], φ = nothing, origin = top_center_pt, rotation = -t.rotation * RotZ{T}(π))
#     # normals of the surfaces show inside the volume primitives.
#     e_top, e_bot, in_mantle, out_mantle
# end

# PartialVaryingTube