%% run_TransCD_IE_WG_TriP0P0RT0.m 
% Trans.conv.-diif.eqn.: Implicit Euler, Weak Galerkin (P0,P0;RT0) 
% James Liu, ColoState; 2012/07--2017/02 

% clc;  
clear all;  close all;
format compact;  format long e;

%% Setting up the problem 
xa = 0;  xb = 1;  yc = 0;  yd = 1;
BndryDescMat = [xa,yc,xb,yc, 0,-1;...
                xb,yc,xb,yd, 1, 0;... 
                xb,yd,xa,yd, 0, 1;...
                xa,yd,xa,yc,-1, 0];
% "0" for convenience 
EqnBIC = EqnBIC_TransCD_Ex11MovingGauss;
% EqnBIC = EqnBIC_TransCD_Ex12MovingGauss;
% EqnBIC = EqnBIC_TransCD_Ex21RotGauss;
% EqnBIC = EqnBIC_TransCD_Ex22RotGauss;
% EqnBIC = EqnBIC_TransCD_Ex23RotCylin;
BndryCondType = [0;1;1;1;1];

%% Discretization: Mesh generation and preparation 
% disp('Generating & preparing a uniform triangular mesh ...');
n = 16 
nx = n;  ny = n;
TriMesh = RectDom_TriMesh_GenUnfm(xa,xb,nx,yc,yd,ny,2);
TriMesh = TriMesh_Enrich1(TriMesh);
TriMesh = TriMesh_Enrich2(TriMesh,BndryDescMat);
TriMesh = TriMesh_Enrich3(TriMesh,BndryDescMat);

%% Setting up Gaussian quadratures 
GAUSSQUAD = SetGaussQuad(5,25,13);

%% Sorting out boundary edges: Dirichlet, Neumann 
DirichletEdge = find(BndryCondType(TriMesh.BndryEdge+1)==1);
NeumannEdge   = find(BndryCondType(TriMesh.BndryEdge+1)==2);

%% Timing 
T = EqnBIC.T;  
NT = 2*n;  MT = NT;  Deltat = T/NT;  tm = (0:NT)*Deltat;
tyme.T = T;  tyme.NT = NT;  tyme.MT = MT; 
tyme.Deltat = Deltat;  tyme.tm = tm;

%% Assembling and solving 
disp('IE + WG(P0,P0;RT0) solving...');
sln = TransCD_IE_WG_TriP0P0RT0_AsmSlv(...
  EqnBIC,TriMesh,tyme,DirichletEdge,NeumannEdge,GAUSSQUAD);

%% Retrieving 
NumerCntrnT = sln(1:TriMesh.NumEms,MT+1);

%% Projecting the initial concentration 
ProjCntrn0 = WG_TriP0P0RT0_ProjScaFxn0(TriMesh,GAUSSQUAD,EqnBIC.fxnc0);

%% Projecting the exact concentration at final time (assuming known) 
ProjCntrnT = WG_TriP0P0RT0_ProjScaFxn1(...
  TriMesh,GAUSSQUAD,EqnBIC.fxnc,(T*MT/NT));

%% Presenting numerical and graphical results 
[min(NumerCntrnT), max(NumerCntrnT)]
% show_TriMesh_ScaEm_ClrImg(TriMesh, ProjCntrn0.em, ...
%   11, 'Projected initial concentration');
show_TriMesh_ScaEm_ClrImg(TriMesh, ProjCntrnT.em, ...
  12, 'Projected concentration at final time');
show_TriMesh_ScaEm_ClrImg(TriMesh, NumerCntrnT, ...
  21, 'Numerical concentration at final time');

return;