%% run_Darcy_WG_TriP0P1BDM1.m 
% Assuming PermKs is an elementwise constant scalar 
% James Liu, ColoState; 2012/07--2018/06 

% clc;  close all;
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_Poisson_WihlerRivieJSC2011_Ex;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Darcy_DurlofWRR1994;
% 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_Enrich1(TriMesh);
TriMesh = TriMesh_Enrich2(TriMesh,BndryDescMat);
TriMesh = TriMesh_Enrich3(TriMesh,BndryDescMat);

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

%% Sampling the permeability 
PermKs = Darcy_SmplnPermSca_TriMesh(EqnBC.fxnKs, TriMesh, GAUSSQUAD);

%% Assembling and Solving: WG(P0,P1;BDM1) 
% tic;
[sln,GlbMat,GlbRHS] = Darcy_WG_TriP0P1BDM1s_AsmSlv(...
  EqnBC,BndryCondType,TriMesh,PermKs,GAUSSQUAD);
% toc

%% Computing numerical  
[NumerPresEm,NumerVelCofBDM1] = Darcy_WG_TriP0P1BDM1s_PresVelFlux(...
  TriMesh,PermKs,sln,GAUSSQUAD);
CofRT0 = TriMesh_ConvBDM1toRT0(TriMesh,NumerVelCofBDM1);
NumerVelEmCntr = TriMesh_ConvRT0toRT00(TriMesh,CofRT0);

%% Computing errors when an exact solution is known 
[L2ErrPres,L2ErrVel] = Darcy_WG_TriP0P1BDM1s_Err(...
  EqnBC,TriMesh,PermKs,sln,GAUSSQUAD);

%% Projecting exact pressure and velocity 
[ProjPresEm,ProjPresEgCof,ProjVelCofBDM1] = ...
  Darcy_WG_TriP0P1BDM1_ProjPresVel(EqnBC,TriMesh,GAUSSQUAD);

%% Presenting
[min(NumerPresEm),max(NumerPresEm)];
[L2ErrPres,L2ErrVel]
show_TriMesh_ScaVecEm_mix(TriMesh, NumerPresEm, NumerVelEmCntr, ... 
  23, 'WG: Numerical pressure and velocity', 1);

return;