%% run_Darcy_WG_PlygnP0P0CW0.m
% Solving Darcy eqn. by WG(P0,P0;CW0) on a convex polygonal mesh 
% James Liu, Zhuoran Wang, ColoState; 2012/07--2018/12 

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

%% Setting up the problem 
% A domain & its boundary (start/ending points, outward unit normal vecs.) 
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];
% Boundary condition type: 1 Dirichlet; 2 Neumann; 0 convenience 
% EqnBC = EqnBC_Poisson_Ex11sinsin;
% BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Poisson_ArnoldBoffiFalkMCOM2002_Ex1;
% BndryCondType = [0;1;1;1;1];
EqnBC = EqnBC_Poisson_MuWangYeIJNAM2015_Case2;
BndryCondType = [0;1;2;1;1];
% EqnBC = EqnBC_Poisson_WheeYotovSINUM2006_Ex2;
% 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 
% % Using PolyMesher 
% tic;
% n = 3
% [Node,Element,Supp,Load,P] = PolyMesher(@UnitSquare, (2^n)^2, n*500);
% PlygnMesh = PlygnMesh_ConvFrom_PolyMesher(Node,Element);
% toc
% 
% Generating a trapezoidal (quadrilateral) mesh 
n = 2^4;  nx = n;  ny= n;  alpha = 1/3;
QuadriMesh = RectDom_QuadriMesh_GenTrapezoid(xa,xb,nx,yc,yd,ny,alpha,2);
PlygnMesh = PlygnMesh_ConvFrom_QuadriMesh(QuadriMesh);
%
% % For a randomly perturbed rough quadrilateral mesh 
% n = 2^3;  nx = n;  ny = n;  delta = 0.25;  status = 2;
% QuadriMesh = RectDom_QuadriMesh_GenLogiRectPtbRand(...
%   xa,xb,nx,yc,yd,ny,delta,status);
% PlygnMesh = PlygnMesh_ConvFrom_QuadriMesh(QuadriMesh);
% 
% % For a rectangular or smooth quadrilateral mesh 
% n = 20;  nx = n;  ny = n;  delta = 0;  status = 2;
% QuadriMesh = RectDom_QuadriMesh_GenLogiRectPtbSmth(...
%   xa,xb,nx,yc,yd,ny,delta,status);
% PlygnMesh = PlygnMesh_ConvFrom_QuadriMesh(QuadriMesh);
% 
PlygnMesh = PlygnMesh_Enrich1(PlygnMesh);
PlygnMesh = PlygnMesh_Enrich2(PlygnMesh,BndryDescMat);
diam = max(PlygnMesh.radius) + min(PlygnMesh.radius);
% [PlygnMesh.NumEms, diam];
% PlygnMesh.maxn 
if PlygnMesh.NumEms<=2^10
  show_PlygnMesh_labels(PlygnMesh,11,0,0,0);
end

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

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

%% WG(P0,P0;CW0)Plygn: Assembling and Solving 
% tic;
sln = Darcy_WG_PlygnP0P0CW0_AsmSlv( ...
  EqnBC, BndryCondType, PlygnMesh, PermK, GAUSSQUAD);
% toc

%% Computing numerical pressure, velocity, flux, ... 
[NumerPresCofP0, NumerVelCofCW0, NumerFlux] = ...
  Darcy_WG_PlygnP0P0CW0_PresVelFlux(PlygnMesh,PermK,sln,GAUSSQUAD);

%% For presentation 
[NumerPresEmAvg, NumerVelEmCntr] = ...
  Darcy_WG_PlygnP0P0CW0_Presn(PlygnMesh, NumerPresCofP0, NumerVelCofCW0);

%% Computing errors when an exact solution is known 
[L2ErrPres,L2ErrVel,L2ErrFlux,L2ErrDiv] = Darcy_WG_PlygnP0P0CW0_Err( ...
  EqnBC, PlygnMesh, NumerPresCofP0, NumerVelCofCW0, GAUSSQUAD);

%% Presenting numerical and graphical results
% [L2ErrPres, min(NumerPresEm), max(NumerPresEm)]
[L2ErrPres, L2ErrVel, L2ErrFlux, L2ErrDiv]
show_PlygnMesh_ScaVecEm_mix(PlygnMesh,NumerPresEmAvg,NumerVelEmCntr, ...
  23,'WG(P_0,P_0;CW_0) numerical pressure & velocity',1);  % Magnified 1X

return;