%% run_DarcyGen_WG_QuadriP0P0AC0.m
% Darcy eqn. in a general form solved by WG(P0,P0;AC0);
% Note: AC0 = RT[0] on rectangles 
% James Liu, ColoState; 2017/07--2021/02 

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

%% Setting up the problem to be solved 
% Boundary condition types: 1 Dirichlet, 2 Neumann, 0 convenience 
EqnBC = EqnBC_DarcyGen_GiraultVassiYotovNumerMath2014_Ex;
BndryCondType = [0; 1; 1; 2; 2];
% The domain & its boundary (start/ending points, outward unit normals) 
xa = 0;  xb = 1;  yc = 0;  yd = 0.5;
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];

%% Mesh generation and preparation 
% % For a rectangular or smooth quadrilateral mesh 
% n = 2^6;  nx = n;  ny= n;  delta = 0.1;  status = 2;
% QuadriMesh = RectDom_QuadriMesh_GenLogiRectPtbSmth(...
%   xa,xb,nx,yc,yd,ny,delta,status);
% For a randomly perturbed rough quadrilateral mesh 
n = 2^6;  nx = n;  ny = n;  delta = 0.25;  status = 2;
QuadriMesh = RectDom_QuadriMesh_GenLogiRectPtbRand(...
  xa,xb,nx,yc,yd,ny,delta,status);
% % For a trapezoidal mesh in ArnoldBoffiFalk_MCOM_2002 
% n = 2^6;  nx = n;  ny= n;  alpha = 0.25;
% QuadriMesh = RectDom_QuadriMesh_GenTrapezoid(xa,xb,nx,yc,yd,ny,alpha,2);
QuadriMesh = QuadriMesh_Enrich2(QuadriMesh,BndryDescMat);
QuadriMesh = QuadriMesh_Enrich3(QuadriMesh,BndryDescMat);
QuadriMesh = QuadriMesh_Enrich4(QuadriMesh,BndryDescMat);
if QuadriMesh.NumEms<=2^6
  show_QuadriMesh_labels(QuadriMesh,11,1,1,1);
end

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

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

%% DarcyGen solved by WG(P0,P0;AC0) 
[sln,CofQhg] = DarcyGen_WG_QuadriP0P0AC0_AsmSlv( ... 
  EqnBC, BndryCondType, QuadriMesh, PermK, GAUSSQUAD);

%% Post-processing: Computing quantities for later use  
[NumerPresCofP0P0, NumerVelCofAC0] ... 
  = DarcyGen_WG_QuadriP0P0AC0_PresVelFlux( ...
    EqnBC, QuadriMesh, PermK, sln, CofQhg, GAUSSQUAD);

%% Further processing for presentation 
[NumerPresEmAvg, NumerVelEmCntr] ... 
  = DarcyGen_WG_QuadriP0P0AC0_Presn( ... 
    QuadriMesh, NumerPresCofP0P0, NumerVelCofAC0);

%% Computing errors when an exact solution is known 
[L2ErrPres,L2ErrVel,L2ErrFlux,L2ErrDiv] ... 
  = DarcyGen_WG_QuadriP0P0AC0_Err( ...
    EqnBC, QuadriMesh, NumerPresCofP0P0, NumerVelCofAC0, GAUSSQUAD);

%% Presenting numerical and graphical results
[L2ErrPres, L2ErrVel]
show_QuadriMesh_ScaVecEm_mix(QuadriMesh, NumerPresEmAvg, NumerVelEmCntr,...
  23, 'DarcyGen: WG numer. pressure & velocity', 1);

return;