%% run_Darcy_MFEM_TriBDM1P0.m 
% Solving Darcy eqn. by MFEM (BDM1,P0) on a triangular mesh 
% Assuming PermKs is an elementwise constant scalar 
% James Liu, ColoState; 2012/07--2018/06 

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

%% Setting up the 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];
EqnBC = EqnBC_Poisson_Ex11sinsin;
BndryCondType = [0;1;1;1;1];
% EqnBC = EqnBC_Darcy_Ex13HetePerm;
% BndryCondType = [0;2;1;2;1];  

%% Generating & preparing a triangular mesh 
n = 2^6 
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 quadratures 
GAUSSQUAD = SetGaussQuad(5,25,13);

%% Sampling the permeability (as a scalar) 
PermKs = Darcy_SmplnPermSca_TriMesh(EqnBC.fxnKs, TriMesh, GAUSSQUAD);

%% Assembling and Solving... 
% tic;
[sln,MatA,MatB,GlbRHS] = Darcy_MFEM_TriBDM1P0s_AsmSlv(...
  EqnBC,BndryCondType,TriMesh,PermKs,GAUSSQUAD);
% toc

%% Computing numerical pressure, velocity, etc. 
[NumerPresEm,CofBDM1NumerVel] = Darcy_MFEM_TriBDM1P0_PresVelFlux(TriMesh,sln);
CofRT0 = TriMesh_ConvBDM1toRT0(TriMesh,CofBDM1NumerVel);
NumerVelEmCntr = TriMesh_ConvRT0toRT00(TriMesh,CofRT0);

%% Computing errors when an exact solution is known 
[L2ErrPres,L2ErrVel] = Darcy_MFEM_TriBDM1P0_Err(EqnBC,TriMesh,sln,GAUSSQUAD);

%% Presenting numerical and graphical results
[L2ErrPres, L2ErrVel, min(NumerPresEm), max(NumerPresEm)]
show_TriMesh_ScaVecEm_mix(TriMesh, NumerPresEm, NumerVelEmCntr,...
  23, 'MFEM: Numerical pressure and velocity', 1);

return;