function [NumerPresEm,NumerVelCofBDM1,NumerFlux] = ...
  Darcy_MFEM_TriBDM1P0_PresVelFlux(TriMesh,sln)
%% Darcy: MFEM(BDM1,P0): Computing numerical pressure, velocity, flux 
% Assuming PermKs is an elementwise constant scalar 
% James Liu, ColoState; 2012/07--2017/02 

%% Mesh info 
NumEms = TriMesh.NumEms;
NumEgs = TriMesh.NumEgs;

%% "Computing" numerical pressure 
NumerPresEm = sln(2*NumEgs+(1:NumEms));

%% Computing numerical velocity 
NumerVelCofBDM1 = zeros(NumEms,6);
for j=1:3
  NumerVelCofBDM1(:,2*j-1) = sln(2*TriMesh.elem2edge(:,j)-1)...
    .*TriMesh.SignEmEg(:,j);
  NumerVelCofBDM1(:,2*j)   = sln(2*TriMesh.elem2edge(:,j))...
    .*TriMesh.SignEmEg(:,j);
end
% Adjusting: swapping 
for ie=1:NumEms
  for j=1:3
    if (TriMesh.SignEmEg(ie,j)==-1) 
      tmp = NumerVelCofBDM1(ie,2*j-1);
      NumerVelCofBDM1(ie,2*j-1) = NumerVelCofBDM1(ie,2*j);
      NumerVelCofBDM1(ie,2*j) = tmp;
    end
  end
end

%% JL20170212: TO BE FINISHED/REVISED 
%% Computing numerical flux 
NumerFlux = zeros(NumEms,3);

return;