% Given: Node,Edge,Face
% asembling of all element stiffness matrices (per element),
% the global stiffness matrix (Neumann everywhere)
% K(nnode,nnode) and the right-hand-side vector

function [K,Face_estiff,Face_eorder] = ...
      Assembly_fine_trianCDR(Node,Edge,Face,Face_estiff,Face_eorder,...
		           Node_flagx,Node_flagy,bndry_edge,...
                           Face_flag,Face_thick,Disco,Discoef,...
                           wh,ndof,lvl_total)

global epsilon flag tita c

nface = size(Face,2);   % number of finite elements
nedge = size(Edge,2);   % number of edges
nnode = size(Node,2);   % number of nodes
K   = spalloc(nnode,nnode,ndof*2*nnode);
rhs= zeros(2*nnode,1);
nall = 2*nnode;
bndry_edge_normal=spalloc(2,nedge,2*nedge);

disp('Begin allocating memory...')
lengthK=nface*ndof*ndof;
KI  = zeros(lengthK,1);
KJ  = zeros(lengthK,1);
KV  = zeros(lengthK,1);
nextK = 0;
disp('...end allocating memory.')

[Gauss_point,Gauss_weight] = Integr_weights_trian;
DER    = shape_der_trian;                    % DER(2x3) %linear b.f.
for iface=1:nface,

   local_node_list=order_face_nodes_tria1(Node,Edge,Face,iface);

   Coord(1:ndof,1) = Node(1,local_node_list(1,1:ndof))';  % Coord(ndof,2)
   Coord(1:ndof,2) = Node(2,local_node_list(1,1:ndof))';

% ------------------------ determine nju 
nju = Discoef(1,Face_flag(iface,1)); % the coefficient jump

% - - - - - generation of the element matrices:
[L_elem,C_elem,M_elem] = Assm_CDR_trian(Gauss_point,Gauss_weight,Coord,DER,wh);
%L_elem = nju*L_elem; % discontinuous coefficients
E_elem = nju*L_elem+C_elem+c*M_elem;

%off Face_eorder(1:ndof,iface,lvl_total)        = local_node_list';
%off Face_estiff(1:ndof,1:ndof,iface,lvl_total) = E_elem;

    for i=1:ndof,
        is = local_node_list(i);
       for j=1:ndof,
          js = local_node_list(j);
%          K(is,js) = K(is,js) + E_elem(i,j);          
	  nextK = nextK + 1;
	  KI(nextK) = is;
	  KJ(nextK) = js;
	  KV(nextK) = E_elem(i,j);
       end
%       rhs(is,1) = rhs(is,1) + rhs_elem(i,1);
    end
    
end             % end iface-loop
%disp(['ququ ' num2str(toc)])

K  = sparse(KI, KJ, KV); %K = (K+K')*0.5;

%disp(['nextK'  int2str(nextK)]) 
%disp(['nextC0' int2str(nextC0)])

% Impose Dirichlet b.c. on the fime mesh matrix

[K,rhs]=Dirichlet(K,rhs,Node,Node_flagx,Node_flagy,nnode);

return