%% run_Darcy_CG_TriP1.m 
% Darcy: Continuous Galerkin (CG) P1 FEM on a triangular 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;

%% Setting up a Darcy flow problem
% A rectangular domain: 4 corners, 4 sides w/ outward unit normal vectors 
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 = 20  % For EqnBC_Darcy_Ex13HetePerm, n needs to be a multiple of 20 
nx = n;  ny = n;
TriMesh = RectDom_TriMesh_GenUnfm(xa,xb,nx,yc,yd,ny,2);
TriMesh = TriMesh_Enrich2(TriMesh,BndryDescMat);
TriMesh = TriMesh_Enrich3(TriMesh,BndryDescMat);
if TriMesh.NumEms<=100 
  show_TriMesh_labels(TriMesh,11,1,1,1);
end

%% Setting up Gaussian quadratures on edges, rectangles, triangles 
GAUSSQUAD = SetGaussQuad(5,25,9);

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

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

%% Assembling and Solving: CG_TriP1 
% tic;
[sln] = Darcy_CG_TriP1_AsmSlv(... 
  EqnBC, TriMesh, PermK, DirichletEdge, NeumannEdge, GAUSSQUAD);
% toc

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

%% Presenting numerical and graphical results 
show_TriMesh_ScaNd_trisurf(TriMesh,sln, ... 
  21, 'CG.TriP1: Numerical pressure');
show_TriMesh_ScaVecEm_mix(TriMesh, NumerPresEm, NumerVelEmCntr, ...
  23, 'CG.TriP1: Numerical pressure and velocity', 1);
show_scalar(FluxDscp, ... 
  24, 'CG.TriP1: Flux discrepancy across edges');
% show_TriMesh_ScaEm_ClrImg(TriMesh, LMCR, ...
%   25, 'CG.TriP1: Local-mass-conservation residue');

return;