%% run_Darcy_FVMD_TriP1.m 
% Solving Darcy eqn. by finite volume method (node-based): tri.mesh, P1 bas.fxn. 
% Assuming PermK is an elementwise constant scalar 
% James Liu, ColoState; 2012/07--2018/06

% clc; 
clear all;  close all; 
format compact;  format long 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];

%% Discretization: Mesh generation and preparation 
n = 2^5 
nx = n;  ny = n;
TriMesh = RectDom_TriMesh_GenUnfm(xa,xb,nx,yc,yd,ny,2);
TriMesh = TriMesh_Enrich1(TriMesh);
TriMesh = TriMesh_Enrich2(TriMesh,BndryDescMat);
TriMesh = TriMesh_Enrich3(TriMesh,BndryDescMat);

%% Setting up quadratures 
GAUSSQUAD = SetGaussQuad(5,25,13);

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

%% Assembling and solving... 
% tic;
[sln,GlbMat,GlbRHS] = Darcy_FVMD_TriP1s_AsmSlv(...
  EqnBC,BndryCondType,TriMesh,PermKs,GAUSSQUAD);
% toc

%% Presenting numerical & graphical results 
% [min(sln), max(sln)]
1-max(sln)
show_TriMesh_ScaNd_trisurf(TriMesh, sln, ...
  21, 'Numerical pressure');

return;