%% run_Darcy_CG_RectQ1.m 
% Darcy: Continuous Galerkin (CG) Q1 FEM on a rectangular mesh 
% Assuming PermK (permeability) is an elementwise constant 2x2 SPD matrix 
% James Liu, ColoState; 2012/07--2018/12 

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

%% Setting up a Darcy flow problem
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];
% 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  % For EqnBC_Darcy_Ex13HetePerm, n needs to be a multiple of 20 
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 tensor (2x2 SPD matrix)
PermK = Darcy_SmplnPerm_RectMesh(EqnBC.fxnK, RectMesh, GAUSSQUAD);

%% Sorting out boundary edges: 1 for Dirichlet, 2 for Neumann 
DirichletEdge = find(BndryCondType(RectMesh.BndryEdge+1)==1);
NeumannEdge   = find(BndryCondType(RectMesh.BndryEdge+1)==2);

%% Assembling and Solving: CG_RectQ1 
% tic;
[sln] = Darcy_CG_RectQ1_AsmSlv(... 
  EqnBC, RectMesh, PermK, DirichletEdge, NeumannEdge, GAUSSQUAD);
% toc

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

%% Computing errors against the known exact solution
err = Darcy_CG_RectQ1_Err(EqnBC,RectMesh,PermK,sln,GAUSSQUAD);

%% Presenting numeircal and graphical results 
err(5)
show_RectMesh_ScaVecEm_mix(RectMesh, NumerPresEm, NumerVelEmCntr, ...
  23, 'CG.RectQ1: Numerical pressure and velocity', 1); 
show_scalar(FluxDscp, ... 
  24, 'CG.RectQ1: Flux discrepancy across edges');

return;