%% run_Darcy_WG_RectQ0Q0RT0.m
% Darcy: WG(Q0,Q0;RT[0]) on a rectangular mesh
% Assuming PermK is an elementwise constant 2x2 SPD matrix 
% James Liu, Farrah Sadre-Marandi, Zhuoran Wang, 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];
% Setting up the Darcy equation & boundary conditions on four sides: 
% 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];
% EqnBC = EqnBC_Darcy_StanfordPerm;
% 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;
RectMesh = RectDom_RectMesh_GenUnfm(xa,xb,nx,yc,yd,ny,2);
RectMesh = RectMesh_Enrich2(RectMesh,BndryDescMat);
RectMesh = RectMesh_Enrich3(RectMesh,BndryDescMat);

%% Setting up Gaussian quadratures on edges, rectangles, triangles 
GAUSSQUAD = SetGaussQuad(5,25,9);

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

%% WG: RectQ0Q0RT0: Assembling and Solving 
% tic;
[sln] = Darcy_WG_RectQ0Q0RT0_AsmSlv( ...
  EqnBC, BndryCondType, RectMesh, PermK, GAUSSQUAD);
% toc

%% Computing numerical pressure, velocity, flux, 
% flux discrepancy, local-mass-conservation residual (LMCR) 
[NumerPresEm, NumerVelCofRT0, NumerFlux, FluxDscp, LMCR] = ... 
  Darcy_WG_RectQ0Q0RT0_PresVelFlux(RectMesh, PermK, sln);
NumerVelEmCntr = NumerVelCofRT0(:,1:2);

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

%% Presenting numerical and graphical results
[min(NumerPresEm), max(NumerPresEm), L2ErrPres, L2ErrVel] 
show_RectMesh_ScaVecEm_mix(RectMesh, NumerPresEm, NumerVelEmCntr, ...
  23, 'WG(Q0,Q0;RT[0]): Numer. pressure & velocity', 2);  % Magnified 2X 
% print -f23 -depsc2 Darcy_Ex13HetePerm_WG_RectQ0Q0RT0_n.eps 
show_scalar(FluxDscp, ... 
  24, 'WG(Q0,Q0;RT[0]): Numer. flux discrepancy across edges');
% show_RectMesh_ScaEm_ClrImg(RectMesh, LMCR, ...
%   25, 'WG(Q0,Q0;RT[0]): Local-mass-conservation residue');

% %% Presenting extra/special results 
% figure(31);
% x = linspace(xa, xb, n);
% y = linspace(yc, yd, n);
% cv = [0.0 : 0.1 : 1.50];
% contour(x, y, reshape(NumerPresEm,n,n), cv);
% xlim([xa xb]);  ylim([yc yd]);
% set(gca, 'fontsize', 20);
% title('WG(Q0,Q0;RT[0]): Numerical pressure contours', 'fontsize', 20);
% axis square;
% % print -f31 -deps2c Ex_PresCntur_WGQ0Q0RT0_n.eps;

return;