%% run_SteadyCD_CG_TriP1.m 
% Solving a steady-state convection-diffusion eqn. by CG P1 on triangular mesh 
% 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 
EqnBC = EqnBC_SteadyCD_Ex21ConstVelBndryLayers;
BndryCondType = [0;1;1;1;1];

%% Discretization: Generating a uniform triangular mesh
n = 32;
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 for SteadyCD: Dirichlet, Neumann 
DirichletEdge = find(BndryCondType(TriMesh.BndryEdge+1)==1);
NeumannEdge   = find(BndryCondType(TriMesh.BndryEdge+1)==2);

%% JL2015111: SPECIAL CASE v=(1,0)
CofRT0NmlzBas = zeros(TriMesh.NumEms,3);
CofRT0NmlzBas(:,1) = 1;  
CofRT0NmlzBas(:,2) = 0;
CofRT0EgBas = TriMesh_RT0_NmlzBas_EgBas(TriMesh,CofRT0NmlzBas);
NmlFlux = zeros(TriMesh.NumEms,3);

%% Solving 
sln = SteadyCD_CG_TriP1_AsmSlv(EqnBC,TriMesh,CofRT0EgBas,NmlFlux,...
  DirichletEdge,NeumannEdge,GAUSSQUAD);

%% Presenting
[min(sln),max(sln)]
show_TriMesh_ScaNd_trisurf(TriMesh, sln, ...
  21, 'Numerical concentration');

return;