function QuadriMesh = QuadriMesh_Enrich1(QuadriMesh)
%% (Enrich1) Enriching a quadrilateral mesh to ...
% James Liu, Zhuoran Wang, ColoState; 2012/07--2017/02 

%% Coefficients for the bilinear mapping 
k1 = QuadriMesh.elem(:,1);  k2 = QuadriMesh.elem(:,2);  
k3 = QuadriMesh.elem(:,3);  k4 = QuadriMesh.elem(:,4);
x1 = QuadriMesh.node(k1,1);  y1 = QuadriMesh.node(k1,2);
x2 = QuadriMesh.node(k2,1);  y2 = QuadriMesh.node(k2,2);
x3 = QuadriMesh.node(k3,1);  y3 = QuadriMesh.node(k3,2);
x4 = QuadriMesh.node(k4,1);  y4 = QuadriMesh.node(k4,2);
QuadriMesh.area = 0.5*( (x2-x1).*(y3-y1) - (x3-x1).*(y2-y1)...
                      + (x3-x1).*(y4-y1) - (x4-x1).*(y3-y1));
QuadriMesh.EmCntr = [0.25*(x1+x2+x3+x4), 0.25*(y1+y2+y3+y4)];
QuadriMesh.CofA = zeros(QuadriMesh.NumEms,4);       
QuadriMesh.CofA(:,1) = x1;
QuadriMesh.CofA(:,2) = x2-x1;
QuadriMesh.CofA(:,3) = x4-x1;
QuadriMesh.CofA(:,4) = (x1+x3)-(x2+x4);
QuadriMesh.CofB = zeros(QuadriMesh.NumEms,4);
QuadriMesh.CofB(:,1) = y1;
QuadriMesh.CofB(:,2) = y2-y1;
QuadriMesh.CofB(:,3) = y4-y1;
QuadriMesh.CofB(:,4) = (y1+y3)-(y2+y4);

return;