% -------------------------------------------------------------------------
% Generate finite element matrices 
% for the convection-diffusion-reaction problem
% d/dx(epsx du/dx) + d/dy(epsy du/dy) + ...
% v1 du/dx + v2 du/dy + c u = f;
% Parameters: 
% epsx and epsy may be constants and used to describe anisotropy
% epsx and epsy may depend on (x,y) and describe discontinuous coefficients
% (v1,v2) describes the velocity field
% c is the reaction term and may simulate the inverse of the time step
%   for time-dependent problems
% -----------------------------------------------------------------------
%
% Works on any given (coarse) mesh, described as {Node, Edge, Face}   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global epsilon flag tita c


figure('renderer','painters'); % if OpelGL crashes


ndof= 3;  % number of points in the finite element
scale_flag=0;   % scale_flag=1 == scaling on
% set jump
jmp=1;
% ------------------------------ load the very coarse mesh {Node,Edge,Face}
[Node,Edge,Face,Node_flagx,Node_flagy,Edge_flagx,Edge_flagy,...
         Face_flag,Face_thick,Face_Node,...
         Disco,Discoef,wh] = Square_ConvDiff;
%         Disco,Discoef,wh] = Square_disco_hier(jmp);
%         Disco,Discoef,wh] = Square_ConvDiff;

figure(1),clf,Bvisual_mesh(Node,Edge,Face,1,1,1,0,16)
%figure(1),clf,Bvisual_mesh(Node,Edge,Face,1,1,1,3,16)
% -------------------- Input parameters ---------
% set anisotropy
if wh=='jp',
   epsx = 1;epsy = 1; 
else
   epsx = 1e-2; epsy = 1e-2; % define anisotropy
end
epsilon = [epsx 0;0 epsy];

% set the vector fiels
% Possible choices:
% flag=0; - no convection
% flag=1; - a rotating field 
flag = 1;  %input('Vector field flag: ');
if flag==1, tita=2*pi/3; 
else
   tita=1;
end
% set reaction term
c = 0;

levels = input('How many times to refine: ');
%for k=2:levels+1
%    figure(k),clf,Bvisual_mesh(Node,Edge,Face,1,1,1,3,16)
%end

lvl_total = levels + 1;
lvl_coars = 1; %1;
nface_lvl=zeros(lvl_total,1);
nface_lvl(1) = size(Face,2);
nnode_lvl(1) = size(Node,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ------------------------------------ initial refinement 
tic
Face_Parent=[];
Face_Node_cell=cell(lvl_total,1);
Face_Node_cell{1}=Face_Node(:,:,1); clear Face_Node

for lvl=1:levels,
    Face_Node  = zeros(3,nface_lvl(lvl));
[Node,Edge,Face,Node_flagx,Node_flagy,Edge_flagx,Edge_flagy,...
Face_flag,Face_thick,Face_Node,Face_Parent,nface_lvl,nnode_lvl] = ...
     my_Refine_hier(Node,Edge,Face,...
        Node_flagx,Node_flagy,Edge_flagx,Edge_flagy,...
        Face_flag,Face_thick,Face_Node,Face_Parent,...
        nface_lvl,nnode_lvl,lvl);
     nface_lvl(lvl+1)=size(Face,2);
     nnode_lvl(lvl+1)=size(Node,2);        
     Face_Node_cell{lvl+1,1} = Face_Node; clear Face_Node 
    %figure(lvl+1),Bvisual_mesh(Node,Edge,Face,0,0,0,3,16)
end
if size(Node,2)<625,Bvisual_mesh(Node,Edge,Face,0,0,0,1,16),end


%  Create the level permutation vectors

nnode = size(Node,2); % the number of nodes
nedge = size(Edge,2); % the number of edges
nface = size(Face,2); % the number of faces
nall  = nnode;        % number of variables 
%nall  = 2*nnode;      % number of displacement variables 

% ---------------- Create mesh-related matrices -----
Edge_Node = spalloc(nedge,nnode,2*nedge);
Face_Edge = spalloc(nface,nedge,3*nface);
for iedge=1:nedge,
    Edge_Node(iedge,Edge(1,iedge))=1;
    Edge_Node(iedge,Edge(2,iedge))=1;
end
for iface=1:nface,
    Face_Edge(iface,Face(1,iface))=1;
    Face_Edge(iface,Face(2,iface))=1;	 
    Face_Edge(iface,Face(3,iface))=1;
end 

bndry_edge = zeros(nedge,1);
bndry_node = zeros(1,nnode);
bndry_edge = sum(Face_Edge);           % The boundary edges are with sum '1'
for iedge=1:nedge,
    if bndry_edge(iedge)==1,
       bndry_edge(iedge)=find(Face_Edge(:,iedge)); % The face containing iedge
     else
       bndry_edge(iedge)=NaN;
    end
end
    
for iedge=1:nedge, 
    if bndry_edge(iedge)==1,
       bndry_node = bndry_node + Edge_Node(iedge,:);
    end
end

Node_Face    = Edge_Node'*Face_Edge';  % Node-to-Face relation
Node_Face_no = sum(Node_Face')./2;     % to how many faces each node belongs
    
% ---------------- Edge_type -----
Edge_type = zeros(size(Edge,2),1);
disp(['End refinement. Elapsed:   ' num2str(toc)])
disp(['Total number of triangles: ' int2str(nface)])
disp(['Total number of nodes:     ' int2str(nnode)])
disp(' Ready with the mesh.')

%%%%%%%%%%%%%% Matrix generation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Face_estiff=[];Face_eorder=[];

[K,Face_estiff,Face_eorder]=...
     Assembly_fine_trianCDR(Node,Edge,Face,Face_estiff,Face_eorder,...
        Node_flagx, Node_flagy,bndry_edge,...
        Face_flag,Face_thick,Disco,Discoef,...
        wh,ndof,lvl_total);  

sol = cd_sol(Node);  


addpath('/home/maya/student/NLA/AGMG_Yvan_Notay/AGMG_3.0/Matlab')
flag_Laplace=0;
if (c==0)&(flag==0)&(epsx==1)&(epsy==1),  % Laplace
    flag_Laplace=1;
end
if flag_Laplace==1,
    addpath('/home/maya/student/NLA/Applications/Multigrid_Demmel')
    makemgdemo_time
    testfmgv_time

   nun=size(K,1);

   % create a similar solution as for the MG code
   X1=find((Node(1,:)>=0.5)&Node(1,:)<=0.75);
   Y1=find((Node(2,:)>=0.25)&Node(2,:)<=0.5);
   X2=find((Node(1,:)>=0.25)&Node(1,:)<=0.5);
   Y2=find((Node(2,:)>=0.5)&Node(2,:)<=0.75);
   Z1=intersect(X1,Y1);
   Z2=intersect(X2,Y2);
   sol=zeros(nun,1);
   sol(Z1)=-2;
   sol(Z2)=1;
   %plot3(Node(1,:),Node(2,:),sol,'*')
end

rhs=K*sol;
tic, sol_ex=K\rhs;tim_direct=toc;

tic,agmg(K,[],[],[],[],[],[],1);tim_preproc=toc;
tic,[x,fl,relres,it,resvec] = agmg(K,rhs,[],[],[],[],[],2);tim_sol=toc;
tim_tot = tim_preproc+tim_sol; 
disp(['Direct solver time: ' num2str(tim_direct)])
disp(' ')
disp(['AGMG preproc. time: ' num2str(tim_preproc)])
disp(['AGMG solution time: ' num2str(tim_sol)])
disp(['AGMG total time:    ' num2str(tim_tot)])
disp(['AGMG iterations:    ' int2str(it)])

semilogy(resvec,'rd-')
if flag_Laplace==1,
   legend('MG','AMG')
end
title('residual(k)'), 
xlabel('iteration number k'), grid