%% run_TransCDR_ELLAM_CG_RectQ1.m 
% ELLAM: Eulerian-Lagrangian Local Adjointed Method 
% 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];
EqnBIC = EqnBIC_TransCDR_CodinaCMAME1998_Ex;
BndryCondType = [0;1;1;1;1];

%% Discretization: Mesh generation and preparation
% disp('Generating & preparing a uniform rectangular mesh...');
n = 16 
nx = n;  ny = n;
RectMesh = RectDom_RectMesh_GenUnfm(xa,xb,nx,yc,yd,ny,2);
RectMesh = RectMesh_Enrich1(RectMesh);
RectMesh = RectMesh_Enrich2(RectMesh,BndryDescMat);
RectMesh = RectMesh_Enrich3(RectMesh,BndryDescMat);

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

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

%% Timing 
T = EqnBIC.T;  
NT = 4*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;

%% Solving 
disp('ELLAM + CG(RectQ1) solving ...');
tic;
NumerCntrn = TransCDR_ELLAM_CG_RectQ1_AsmSlv(EqnBIC,RectMesh,tyme,...
  DirichletEdge,NeumannEdge,GAUSSQUAD);
toc

%% Presenting numerical and graphical results 
[min(NumerCntrn(:,MT+1)), max(NumerCntrn(:,MT+1))]
% show_RectMesh_ScaNd_mesh(RectMesh, NumerCntrn(:,1), ...
%   21, 'Initial concentration');
show_RectMesh_ScaNd_mesh(RectMesh, NumerCntrn(:,NT+1), ...
  22, 'Final concentration');

return;