%% run_Darcy_MFEM_RectRT0Q0.m 
% Darcy: MFEM(RT[0],Q0) on a rectangular mesh
% Assuming PermKs is an elementwise constant scalar 
% James Liu, ColoState; 2012/07--2018/12 

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

%% Setting up a Darcy flow problem
% Setting up a rect.domain: 4 corners, 4 sides with outward unit normals
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];
% Setting up 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_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;
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 (elementwise constant scalar)
PermKs = Darcy_SmplnPermSca_RectMesh(EqnBC.fxnKs, RectMesh, GAUSSQUAD);

%% MFEM(RT[0],Q0): Assembling and Solving... 
% tic;
[sln] = Darcy_MFEM_RectRT0Q0s_AsmSlv( ...
  EqnBC,BndryCondType,RectMesh,PermKs,GAUSSQUAD);
% toc

%% Computing numerical pressure, velocity, flux, 
% flux discrepancy, local-mass-conservation residual (LMCR) 
[NumerPresEm, NumerVelEmCntr, NumerFlux, FluxDscp, LMCR] = ...
  Darcy_MFEM_RectRT0Q0_PresVelFlux(RectMesh, sln, BndryDescMat);

%% Computing errors assuming an exact solution is known 
[L2ErrPres] = Darcy_MFEM_RectRT0Q0_Err(EqnBC,RectMesh,sln,GAUSSQUAD);

%% Presenting numerical & graphical results 
[min(NumerPresEm), max(NumerPresEm), L2ErrPres] 
show_RectMesh_ScaVecEm_mix(RectMesh, NumerPresEm, NumerVelEmCntr, ...
  23, 'MFEM(RT[0],Q0): Numeri. pressure & velocity', 2);  % Magnified 2X 

return;