function [gMat,gVec,n2d]=FEMEquationM(nodes,elem,edges,aFunc,bFunc,fFunc,gDFunc,gNFunc)
%[...] = FEMEquationM (...)
%  sets up the system of linear equations for a numerical solution of a PDE
%
%  [A,b,n2d]=FEMEquation(nodes,elem,edges,'a','b','f','gD','gN')
%     nodes elem edges  describe the mesh 
%         see ReadMesh() for the description of the format
%     'a','b','f','gD','gN' are the functions and coefficients 
%         for the boundary value problem. They can be given as a scalar value
%         or as a sting with the function name 
%
%  div(a*grad u) - b*u = f     in domain
%                    u = gD    on Dirichlet section of the boundary
%              a*du/dn = gN    on Neumann section of the boundary 
%
%
% A   is the matrix of the system to be solved.
% b   is the RHS of the system to be solved.
% n2d is the renumbering of the nodes to the DOF of the system
%     n2d(k)=0  indicates that node k is a Dirichlet node
%     n2d(k)=nn indicates that the value of the solution at node k 
%               is given by u(nn)
%
%see also ReadMesh, ShowMesh, FEMSolve, ShowSolution, FEMValueM

[n,n2d]=FindDOF(nodes);
gMat=zeros(n,n);
gVec=zeros(n,1);

% insert the element matrices and vectors into the global matrix
for k=1:length(elem)
  corners=[nodes(elem(k,1),1:2);nodes(elem(k,2),1:2);nodes(elem(k,3),1:2)];
  [mat,vec]=ElementContribution(corners,aFunc,bFunc,fFunc);
  dofs=[n2d(elem(k,1)),n2d(elem(k,2)),n2d(elem(k,3))];
  for k1=1:3
  if dofs(k1)>0
    for k2=1:3
      if dofs(k2)>0 
	gMat(dofs(k1),dofs(k2)) = gMat(dofs(k1),dofs(k2)) + mat(k1,k2);
      else
	if is_scalar(gDFunc) gD=gDFunc;
	            else     gD=feval(gDFunc,corners(k2,:));endif
	gVec(dofs(k1)) = gVec(dofs(k1)) +  mat(k1,k2)*gD;
      endif % dofs(k2)>0
    endfor  % k2
    gVec(dofs(k1)) = gVec(dofs(k1)) + vec(k1);
  endif  % dofs(k1)
  endfor % k1
endfor

% insert the edge contributions
for k=1:length(edges)
if edges(k,3)>1
  corners=[nodes(edges(k,1),1:2);nodes(edges(k,2),1:2)];
  edgeVec=ElementContributionEdge(corners,gNFunc);
  dofs=[n2d(edges(k,1)),n2d(edges(k,2))];
   if dofs(1)>0  gVec(dofs(1)) = gVec(dofs(1)) - edgeVec(1); endif
   if dofs(2)>0  gVec(dofs(2)) = gVec(dofs(2)) - edgeVec(2); endif
endif
endfor   % k
endfunction
