%% run_Darcy_MFEM_TriRT0P0.m 
% Solving Darcy eqn. by MFEM(RT0,P0) on a triangular mesh 
% Assuming PermK is an elementwise constant 2*2 SPD matrix 
% James Liu, ColoState; 2012/07--2021/01 

% 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 = 2^4;  % n needs to be a multiple of 20 for EqnBC_Darcy_Ex13HetePerm 
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);

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

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

%% Assembling and solving... 
[sln] = Darcy_MFEM_TriRT0P0_AsmSlv( ...
  EqnBC, BndryCondType, TriMesh, PermK, GAUSSQUAD);

%% Computing numerical pressure, velocity, flux, etc. 
[NumerPresEm, NumerVelEmCntr, NumerFlux, FluxDscp, LMCR] = ...
  Darcy_MFEM_TriRT0P0_PresVelFlux(BndryDescMat, TriMesh, sln);

%% Computing errors 
[L2ErrPres] = Darcy_MFEM_TriRT0P0_Err(EqnBC,TriMesh,sln,GAUSSQUAD);

%% Presenting numerical & graphical results 
[min(NumerPresEm), max(NumerPresEm), L2ErrPres]
show_TriMesh_ScaVecEm_mix(TriMesh, NumerPresEm, NumerVelEmCntr, ...
  23, 'MFEM(RT0,P0): Numeri. pressure & velocity', 1);  % Magnified 1X 

return;