### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ f6d978ea-b349-11eb-309e-f5a517efc907
begin 
	using LinearAlgebra 
	using Plots 
	using QuadGK
end

# ╔═╡ 99268461-9e3f-4685-b530-e24febfb539e
#= 
IDEA OF MULTIGRID : 
		
	G is a mesh of a 1D domain. 
	G1 is another, but coarser. 
	G1's operator has higher frequency eigenvalues as opposed to  G. 
	Use G1 to weed out these values. 
	Use smoothers alternatively. 
=# 

# ╔═╡ b422b21d-b9d3-4aad-b29b-d0669331abc8
 
begin
	e = 5; #  Tune e for different basis functions 
	N = 10;
	Ni = 100; 
	xmin = 0;
	xmax = 1; 
	x = LinRange(xmin,xmax,N);
	xi = LinRange(x[e-1],x[e+1],Ni);
	h = x[2] - x[1]; 
end

# ╔═╡ b11596d4-ee5d-41ec-af8f-96734ba2f579
# This shows a heirarchy of grids	
grid1 = plot(x,zeros(length(x)),marker="*")

# ╔═╡ c1463b68-12e6-4e8c-9902-c2cc7e6caa91
begin
	# Plot the basis functions 
	y = zeros(length(xi)); 
	for i = 1:length(xi);  
		if(xi[i]>=xi[1] && xi[i]<=x[e])
			y[i] = 1/h * (xi[i]-x[e-1]);	
		elseif (xi[i]>x[e] && xi[i]<x[e+1]) 
			y[i] = 1/h * (x[e+1]-xi[i]);
		end
	end
	plot!(xi,y,title="Linear Basis Functions");
end

# ╔═╡ d5da54e6-0a4f-4f43-a7ab-8791574b32c6
# Now we use the linear basis functions and obtain an approximation of the Laplacian operator. 

begin 
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	function θ(x,i,N)
		n = length(x); 
		y = zeros(n); 
		if i==1 
			for k=1:N 
				if x[k]<= x[i+1] 
					y[k] = (N-1)*(x[i+1]-x[k]);
				end
			end
		elseif i==N 
			for l = 1:N 
				if x[l] >= x[i-1] 
					y[l] = (N-1) * (x[l]-x[i]);
				end
			end
		else 
			for j=2:n-1 
				if(x[j]>=x[i-1] && x[j]<=x[i])
					y[j] = (N-1) * (x[j]-x[i-1]);
				elseif(x[j]>=x[i] && x[j]<=x[i+1])
					y[j] = (N-1) * (x[i+1]-x[j]); 
				end
			end
		end
		return y;
	end
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	function θ_scalar(x,i,xi,N)
			if(x>=xi[i-1] && x<=xi[i])
				y = (N-1) * (x-xi[i-1]);
			elseif(x>=xi[i] && x<=xi[i+1])
				y = (N-1) * (xi[i+1]-x); 
			else
				y=0;
			end
		return y;
	end
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	function stiffness_matrix(n)
		md = 2.0*ones(n-2);
		ld = -1.0*ones(n-3); 
		ud = -1.0*ones(n-3);
		return (n-1)*Tridiagonal(ld,md,ud); 
	end
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	function intergrate(x,i,N) 
		domain = LinRange(0,1,N+2);
		return sin(pi*x) * θ_scalar(x,i,domain,N);  
	end
	
	function compute_rhs(n)
		n = n-2; 
		integral = zeros(n);	
		for i=1:length(integral)  
		integral[i] = 1.5*quadgk(q->intergrate(q,i+1,n),0,1,rtol=1e-8)[1];
		end
		return integral; 
	end
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
end

# ╔═╡ a5f1a911-3f71-4630-a872-99514d52e5f3
begin
	nk = 12;
	a = LinRange(xmin,xmax,nk); 
	z = θ(a,3,nk);
	plot(a,z)
end

# ╔═╡ feba23c9-1b66-400e-9ffb-0dc5fc65896f
begin
	# Get the stiffness matrix and the RHS : The RHS is evaluated and plotted. 
	n = 7;
	test_A = stiffness_matrix(n)
	test_b = compute_rhs(n);
	plot(test_b)
	#test_A
	#test_b
end

# ╔═╡ b8b7b670-6039-4bd9-9f95-8652f6e03528
md"""
	
	#NESTED MULTIGRID

	Start on a coarse grid h_0 

	Iterate using a smoother and remove one range of high frequencies. 

	Interpolate solution to a finer grid - h_1
	
	Start step 1 with a smaller time step h_1. 
"""

# ╔═╡ a217bef6-3c9f-4d12-8b3f-eb937987c709
# Description

md"""

To implement a nested multigrid method. 

We would need access to a smoother. The kind of smoother to choose depends on the nature of the matrix A that we solve. 

For the finite element approximation of the poisson equation we have a matrix A that is not diagonally dominant. 

So solve the problem we need a method for which convergence could be proven. 

In comes the Gauss Seidel method. 

"""

# ╔═╡ f2a7eb59-3617-47fe-be19-f6873d0f3c56
function gauss_seidel_relaxation(A,b,x0,max_iterations)
	x = x0; 
	res = b - A*x; 
	for iter = 1:max_iterations 
		for i = 1:length(x)
			x[i] = (1/A[i,i])*(b[i] - dot(A[i,:],x) + A[i,i]*x[i]);
		end
		res = b - A*x; 
	end 
	return x
end

# ╔═╡ abf72f94-406b-41a6-bbdd-fa3ac2919531
# To perform a nested iteration 
begin 
	n_nodes = 7;
	A = stiffness_matrix(n_nodes);
	b = compute_rhs(n_nodes);
	x0 = zeros(n_nodes-2); 
	
	y_appx = gauss_seidel_relaxation(A,b,x0,10); 
	y_appx = vcat(0,y_appx,0);
	println(y_appx)
	
	x_nodes = LinRange(0,1,100); 
	y_actual = zeros(length(x_nodes)); 
	for i=1:length(x_nodes)
		y_actual[i] = 1/pi^2 * sin(pi*x_nodes[i]);
	end
	plot(x_nodes,y_actual,label="Actual Solution");
	plot!(LinRange(0,1,length(y_appx)),y_appx,label="Approximate solution");
end

# ╔═╡ Cell order:
# ╠═f6d978ea-b349-11eb-309e-f5a517efc907
# ╠═99268461-9e3f-4685-b530-e24febfb539e
# ╠═b422b21d-b9d3-4aad-b29b-d0669331abc8
# ╠═b11596d4-ee5d-41ec-af8f-96734ba2f579
# ╠═c1463b68-12e6-4e8c-9902-c2cc7e6caa91
# ╠═d5da54e6-0a4f-4f43-a7ab-8791574b32c6
# ╠═a5f1a911-3f71-4630-a872-99514d52e5f3
# ╠═feba23c9-1b66-400e-9ffb-0dc5fc65896f
# ╟─b8b7b670-6039-4bd9-9f95-8652f6e03528
# ╟─a217bef6-3c9f-4d12-8b3f-eb937987c709
# ╠═f2a7eb59-3617-47fe-be19-f6873d0f3c56
# ╠═abf72f94-406b-41a6-bbdd-fa3ac2919531
