%% run_SteadyCD_FVMDU_TriP1.m 
% 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];
% EqnBC = EqnBC_SteadyCD_Ex10sinsin;
% EqnBC = EqnBC_SteadyCD_Ex11ConstVelBndryLayer;
EqnBC = EqnBC_SteadyCD_Ex13Gauss;
% EqnBC = EqnBC_SteadyCD_Ex21ConstVelBndryLayers;
% EqnBC = EqnBC_SteadyCD_Ex22ConstVelBndryIntrLayers;
% EqnBC = EqnBC_SteadyCD_Ex23RotVelBndryLayer;
% EqnBC = EqnBC_SteadyCD_Ex24RotVelBndryIntrLayers;
BndryCondType = [0;1;1;1;1];

%% Discretization: Mesh generation and preparation 
n = 2^5 
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 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);

%% Assembling and solving... 
tic;
[sln,GlbMat,GlbRHS] = SteadyCD_FVMDU_TriP1_AsmSlv(EqnBC,...
  TriMesh,DirichletEdge,NeumannEdge,GAUSSQUAD);
toc

%% "Computing" many things
NumerCntrnNd = sln(1:TriMesh.NumNds);
NumerCntrnEm = TriMesh_ConvCP1toP0(TriMesh, NumerCntrnNd);

%% Presenting numerical & graphical results 
[min(sln), max(sln)]
show_TriMesh_ScaNd_trisurf(TriMesh, NumerCntrnNd, ...
  21, 'FVMDU: Numerical concentration');
show_TriMesh_ScaEm_ClrImg(TriMesh, NumerCntrnEm, ...
  22, 'FVMDU: Numerical concentration');

return;