%
%
% 2D streamfunction vorticity Hamiltonian particle FEM code = HPFEM
% By Onno Bokhove April 2010
%
% nabla^2 psi(x,y,t) = omega(x,y,t)
%
%
global xxn yyn Index
%
% Read constants and make mesh: nodes and index file
%
HPFEM10cst;
%
% Assembley of cst S-matrix & bb-vector in: S_ij psi_j = bb_j
%
% loop over all Nk elements
%
S = zeros(ntfa,ntfa); % initialize
RHSv = zeros(ntfa,1); 
psi = zeros(ntfa,1); 
psit = zeros(Nn,1); 
noelements = Nk
%
% Location of main boundary nodes in terms of global node number:
% South boundary: nno = 1:1:Nx
% East boundary: nno = Nx+1:Nx:(Ny-2)*Nx+1
% West boundary: nno = Nx:Nx:(Ny-1)*Nx
% North boundary: nno = (Ny-1)*Nx+1:1:Nx*Ny
%
%
HPFEM10Selem(1,1,1,Nloc,g3,rcoef)
%
%
%
switch nquadratic
case 0 % linear basis functions
 for k = 1:Nk
  % 
  % loop over all local nodes
  %
  for nal = 1:Nloc
    m = Index(nal,k); % global node number of local node nal in element k
    %
    % ONLY DONE FOR LINEAR!
    %
    if ( mnindex(m) > ntfa )
      % m is a Dirichlet boundary node
	???
      %
    else % m is not at Dirichlet boundary node in linear case
      %
      % RHS forcing
      %
	???
      %
      for nbe = 1:Nloc
        l = Index(nbe,k);
        if ( mnindex(l) > ntfa )
           % On boundary;
           % Subtract Matrix values times Dirichlet value from RHsvector
           %
        else % not a boundary node 

        end
      end
    end
  end
 end
case 1 % quadratic basis functions
 for k = 1:Nk
  % 
  % loop over all local nodes
  %
  for nal = 1:Nloc
    m = Index(nal,k); % global node number of local node nal in element k
    %
    if ( mnindex(m) > ntfa )
      % m is a Dirichlet boundary node
      % Set psit vector to Dirichlet value

      %
    else % m interior node; linear case
      %
      % RHS forcing
      %

      %
      for nbe = 1:Nloc
        l = Index(nbe,k);
        if (  mnindex(l) > ntfa )
           % Subtract Matrix values times Dirichlet value from RHsvector

           %
        else % interior node 
           S(mnindex(m),mnindex(l)) = 
        end
      end
    end
  end
 end
end
S  = sparse(S);
assembleyisdone = 1
%
%
% Test numerical solution for psi(x,y,0); i.e., initial condition at t=0
% against exact solution omega(x,y.0) = cos x + cos y
%
%
figure(112); clf
psi = S\RHSv;
psit(1:ntfa) = psi;
xlul = 0:0.01:Lx;
ylul = 0:0.01:Ly;
[XI,YI] = meshgrid(xlul,ylul);
[XI,YI,ZI0] = griddata(xxnm,yynm,psit',XI,YI);
cvec = [0 1 2 4 5 6 7 8 10]*0.0005; % 
cvecm = -[1 2 4 5 6 7 8 10]*0.0005; % 
colorbar;
contourf(XI,YI,ZI0); hold on;
colorbar;
%    
xlabel('x','fontsize',18);    
ylabel('y','fontsize',18); 
%
%
%