%% run_Darcy_WG_TriP0P0RT0.m 
% Darcy: WG(P0,P0;RT0) on a triangular mesh 
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% James Liu, Farrah Sadre-Marandi, ColoState; 2012/07--2018/12 

% clc;  
clear all;  close all;
format compact;  format short e;

%% Setting up a Darcy flow problem
% A rectangular domain: 4 corners, 4 sides w/ outward unit normal vectors 
xa = 0;  xb = 1;  yc = 0;  yd = 1;
% xa = 0;  xb = 5;  yc = 0;  yd = 1;  % For EqnBC_Darcy_StanfordPerm 
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];
% 1 for Dirichlet, 2 for Neumann, 0 for coding convenience in Matlab 
EqnBC = EqnBC_Poisson_Ex11sinsin;
BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_WihlerRivieJSC2011_Ex;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_ArboWheeYotovSINUM1997_Ex2;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Darcy_Ex13HetePerm;
% BndryCondType = [0;2;1;2;1];

%% Discretization: Mesh generation and preparation 
n = 40  % n needs to be a multiple of 20 for EqnBC_Darcy_Ex13HetePerm 
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 Gaussian quadratures on edges, rectangles, triangles 
GAUSSQUAD = SetGaussQuad(5,25,9);

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

%% WG(P0,P0;RT0): Assembling and Solving 
% tic;
[sln] = Darcy_WG_TriP0P0RT0_AsmSlv( ...
  EqnBC, BndryCondType, TriMesh, PermK, GAUSSQUAD);
% toc

%% Computing numerical pressure, velocity, flux, etc. 
[NumerPresEm, NumerVelCofRT0, NumerFlux, FluxDscp, LMCR] = ... 
  Darcy_WG_TriP0P0RT0_PresVelFlux(TriMesh, PermK, sln);
NumerVelEmCntr = NumerVelCofRT0(:,1:2);

%% Computing errors when an exact solution is known 
[L2ErrPres, L2ErrVel] = Darcy_WG_TriP0P0RT0_Err( ...
  EqnBC, TriMesh, NumerPresEm, NumerVelCofRT0, GAUSSQUAD);

%% Presenting numerical & graphical results
[min(NumerPresEm), max(NumerPresEm), L2ErrPres, L2ErrVel]
show_TriMesh_ScaVecEm_mix(TriMesh, NumerPresEm, NumerVelEmCntr, ...
  23, 'WG(P0,P0;RT0): Numerical pressure & velocity', 1);
% show_TriMesh_ScaEm_ClrImg(TriMesh, LMCR, ... 
%   24, 'WG(P0,P0;RT0): Local-mass-conservation residue');

return;