function [elMat,elVec]=ElementContribution(corners,aFunc,bFunc,fFunc)
% [...] =  ElementContribution(...)
%   compute the element stiffness matrix of one element
%
%[elMat,elVec]=ElementContribution(corners,aFunc,bFunc,fFunc)
%  corners  coordinates of the three corners forming the triangular element
%  aFunc bFunc fFunc  function files or scalars for the coefficient functions
%  elMat element stiffness matrix
%  elVec vector contribution to the RHS of the linear equation
 
l1=0.20257301464691267760;  l2=0.94028412821023017954;
w1=0.062969590272413576298; w2=0.066197076394253090369;
w3=0.1125;         w=[w1,w1,w1,w2,w2,w2,w3]';

G = [l1/2 l1/2; 1-l1,l1/2; l1/2, 1-l1;
     l2/2 l2/2; 1-l2,l2/2; l2/2, 1-l2; 1/3 1/3];

T = [corners(2,:)-corners(1,:); corners(3,:)-corners(1,:)];  detT = abs(det(T));
#P=corners(1,:)+G*T
P = G*T; P(:,1) = P(:,1)+corners(1,1); P(:,2) = P(:,2)+corners(1,2);

InterpolationMatrix = [...
 4.743526085855386e-01   -8.076859419188719e-02   -8.076859419188719e-02    4.103582626313830e-02    3.230743767675487e-01     3.230743767675487e-01
   -8.076859419188719e-02    4.743526085855386e-01   -8.076859419188719e-02    3.230743767675487e-01    4.103582626313830e-02     3.230743767675487e-01
   -8.076859419188719e-02   -8.076859419188719e-02    4.743526085855386e-01    3.230743767675487e-01    3.230743767675487e-01     4.103582626313830e-02
   -5.258390110254538e-02   -2.807494322307878e-02   -2.807494322307878e-02    8.841342417640726e-01    1.122997728923152e-01     1.122997728923152e-01
   -2.807494322307878e-02   -5.258390110254534e-02   -2.807494322307878e-02    1.122997728923151e-01    8.841342417640726e-01     1.122997728923151e-01
   -2.807494322307878e-02   -2.807494322307878e-02   -5.258390110254534e-02    1.122997728923151e-01    1.122997728923151e-01     8.841342417640727e-01
   -1.111111111111111e-01   -1.111111111111111e-01   -1.111111111111111e-01    4.444444444444444e-01    4.444444444444445e-01     4.444444444444445e-01];

if ischar(fFunc) val = feval(fFunc,P);  else val = fFunc(:); end%if
elVec = detT*InterpolationMatrix'*(w.*val);

if ischar(bFunc) val = feval(bFunc,P);  else val = bFunc(:); end%if
elMat = detT*InterpolationMatrix'*diag(w.*val)*InterpolationMatrix;

Mxi = [
-2.189707941412349e+00   -5.948539707061746e-01  0  4.051460292938254e-01   -4.051460292938254e-01    2.784561912118524e+00;
+5.948539707061746e-01    2.189707941412349e+00  0  4.051460292938254e-01   -4.051460292938254e-01   -2.784561912118524e+00;
+5.948539707061746e-01   -5.948539707061746e-01  0  3.189707941412349e+00   -3.189707941412349e+00    0.000000000000000e+00;
+7.611365128409209e-01    8.805682564204604e-01  0  1.880568256420460e+00   -1.880568256420460e+00   -1.641704769261381e+00;
-8.805682564204604e-01   -7.611365128409209e-01  0  1.880568256420460e+00   -1.880568256420460e+00    1.641704769261381e+00;
-8.805682564204604e-01    8.805682564204604e-01  0  2.388634871590791e-01   -2.388634871590791e-01    0.000000000000000e+00;
-3.333333333333335e-01    3.333333333333333e-01  0  1.333333333333333e+00   -1.333333333333333e+00    2.220446049250313e-16];

Mnu = [
-2.189707941412349e+00    0  -5.948539707061746e-01    4.051460292938254e-01    2.784561912118524e+00   -4.051460292938254e-01;
 5.948539707061746e-01    0  -5.948539707061746e-01    3.189707941412349e+00    0.000000000000000e+00   -3.189707941412349e+00;
 5.948539707061746e-01    0   2.189707941412349e+00    4.051460292938254e-01   -2.784561912118524e+00   -4.051460292938254e-01;
 7.611365128409209e-01    0   8.805682564204604e-01    1.880568256420460e+00   -1.641704769261381e+00   -1.880568256420460e+00;
-8.805682564204604e-01    0   8.805682564204604e-01    2.388634871590791e-01    0.000000000000000e+00   -2.388634871590791e-01;
-8.805682564204604e-01    0  -7.611365128409209e-01    1.880568256420460e+00    1.641704769261382e+00   -1.880568256420460e+00;
-3.333333333333335e-01    0   3.333333333333333e-01    1.333333333333333e+00    4.440892098500626e-16   -1.333333333333333e+00];

if ischar(aFunc) val = feval(aFunc,P);  else val = aFunc(:); end%if

tt = T(2,2)*Mxi-T(1,2)*Mnu;
Ax = tt'*diag(w.*val)*tt;

tt = -T(2,1)*Mxi+T(1,1)*Mnu;
Ay = tt'*diag(w.*val)*tt;

elMat = elMat + (Ax+Ay)/detT;
end%function
