% Compute the element stiffness matrices corresponding to the
% scalar convection-difuson equation 
% (singularly perturbed, anisotropy and discontinuous coefficients)
% FUN - shape functions for the displacements
% DER - derivatives of FUN
% --------------------------------------------------------------------
function [L_elem,C_elem,M_elem] = Assm_ConvDif_trian(Gauss_point,Gauss_weight,Coord,DER,wh)
     
global epsilon flag tita c
%disp(['Anisotropy: ' num2str(epsilon(1,1)) ',' num2str(epsilon(2,2))])

      nip   = 3;
      nodof = 2;
      ndof  = 6;              % ndof = nip*nodof
      L_elem= zeros(3,3);     % element diffusion matrix
      C_elem= zeros(3,3);     % element convection matrix
      M_elem= zeros(3,3);     % element mass matrix
% ------------------------------------------------------------------
% Gauss_point(2,3)        Coord(3,2)
% FUN(3)                  DER(2,3)          Jac(2,2)     IJac(2,2)
% Deriv(2,3)              
%
      b = b_vecV(Coord,flag,tita);
      for k=1:nip                         % nip = number of integration points
         bb=[b(nip,1) 0;0 b(nip,2)];
         FUN    = shape_fun_trian(Gauss_point,k);     % FUN(1,3)
%off         DER    = shape_der_trian;                    % DER(2x3) %linear b.f.
	     Jac    = DER *Coord;
         Det    = determinant2_m(Jac);
         IJac   = inv(Jac);
         Deriv  = IJac*DER;                         % (2x3)=(2x2)*(2x3)
	     L = Deriv'*(epsilon*Deriv);                  % anisotropic Laplacian
%	     C = b(1)*DER(1,:)'*FUN + b(2)*DER(2,:)'*FUN; % conv-diff
         C = (bb*DER)'*[FUN;FUN];
	     M = FUN'*FUN;                                % mass
         L_elem = L_elem + Det*Gauss_weight(k)*L;
         C_elem = C_elem + Det*Gauss_weight(k)*C;
         M_elem = M_elem + Det*Gauss_weight(k)*M; 
     end
return


%