%% run_Poisson_WG_QuadriP0P0AC0.m
% Solving Poisson eqn. by WG(P0,P0;AC0) on a quadrilateral mesh 
% Assuming CofK is an elementwise constant 2x2 SPD matrix 
% James Liu, ColoState; 2012/07--2021/03 

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

%% Setting up the problem 
% A domain & its boundary (start/ending points, outward unit normal vec.) 
% Boundary condition type: 1 Dirichlet; 2 Neumann; 0 convenience 
EqnBC = EqnBC_Poisson_Ex11sinsin; 
BndryCondType = [0;1;1;1;1];
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];

%% Discretization: Mesh generation and preparation 
% For a rectangular or smooth quadrilateral mesh 
n = 2^5;  nx = n;  ny= n;  delta = 0;  status = 2;
QuadriMesh = RectDom_QuadriMesh_GenLogiRectPtbSmth(...
  xa,xb,nx,yc,yd,ny,delta,status); 
% % 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);
% % For a trapezoidal but asymp. paral. mesh (ArnoldBoffiFalk_MCOM_2002)
% m = 3;  status = 3;
% QuadriMesh = RectDom_QuadriMesh_GenTZasymp(xa,xb,yc,yd,m,status);
% % For a trapezoidal mesh in ArnoldBoffiFalk_MCOM_2002 
% n = 2^6;  nx = n;  ny= n;  alpha = 0.5;
% QuadriMesh = RectDom_QuadriMesh_GenTrapezoid(xa,xb,nx,yc,yd,ny,alpha,2); 
QuadriMesh = QuadriMesh_Enrich2(QuadriMesh,BndryDescMat);
QuadriMesh = QuadriMesh_Enrich4(QuadriMesh,BndryDescMat);
if QuadriMesh.NumEms<=2^8
  show_QuadriMesh_labels(QuadriMesh,11,0,0,0);
end

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

%% Sampling the permeability (2x2 SPD matrix)
CofK = Poisson_SmplnCofK_QuadriMesh(EqnBC.fxnK, QuadriMesh, GAUSSQUAD);

%% WG: QuadriQ0Q0AC0 quadri.mesh: Assembling and Solving 
[sln] = Poisson_WG_QuadriP0P0AC0_AsmSlv(...
  EqnBC,BndryCondType,QuadriMesh,CofK,GAUSSQUAD);

%% Post-processing 
NumEms = QuadriMesh.NumEms;
NumerPresEm = sln(1:NumEms);

%% Presenting numerical and graphical results
[min(NumerPresEm), max(NumerPresEm)]
show_QuadriMesh_ScaEm_ClrImg(QuadriMesh, NumerPresEm, ...
  21, 'WG Numerical pressure');

return;