%% run_Darcy_DG_TriP2.m 
% Assuming PermK is an elementwise scalar 
% James Liu, ColoState; 2012/07--2016/09

% clc;
clear all;  close all;
format compact;  format short 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_Poisson_Ex11sinsin;
BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Darcy_Ex13HetePerm;
% BndryCondType = [0;2;1;2;1];

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

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

%% Sampling the permeability 
PermK = Darcy_SmplnPerm_TriMesh(EqnBC.fxnK,TriMesh,GAUSSQUAD);

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

%% Assembling and solving: DG TriP2 
disp('DG_TriP2: Assembling and solving...');
% Note: beta=-1 for SIPG;  beta=0 for IIPG;  beta=1 for NIPG  
alpha = 10;  beta = 1;
% tic;
[CofNS,GlbMat,GlbRHS] = Darcy_DG_TriP2_AsmSlv(EqnBC,...
  TriMesh,PermK,DirichletEdge,NeumannEdge,GAUSSQUAD,alpha,beta);
% toc

%% Postprocessing 
[PresAvgEm,CofBDM1NumerVel,NumerFlux] = Darcy_DG_TriP2_PostProc(...
  BndryDescMat,TriMesh,PermK,CofNS,GAUSSQUAD);
CofRT0 = TriMesh_ConvBDM1toRT0(TriMesh,CofBDM1NumerVel);
NumerVelEmCntr = TriMesh_ConvRT0toRT00(TriMesh,CofRT0);

% %% Computing errorss 
err = Darcy_DG_TriP2_Err(EqnBC,TriMesh,PermK,CofNS,GAUSSQUAD);
ErrPP = Darcy_DG_TriP2_ErrPP(EqnBC,TriMesh,PermK,...
  CofBDM1NumerVel,GAUSSQUAD);

%% Presenting numerical and graphical results 
[err(2),ErrPP(1)]
% NumerFlux.BndryFlux(2);
show_TriMesh_ScaVecEm_mix(TriMesh, PresAvgEm, NumerVelEmCntr, ...
  23, 'DG: Numerical pressure and velocity', 2);

return;