%
%
% HPFEM10cstst.m : define constants
%
% Use at your own risk; please check and plot the mesh
%
%
nquadratic = 1; % 0: linear basis functions with 3 local nodes per triangle (second order formally);
	    	% 1: quadratic basis functions with 6 local nodes per triangle (third order formally)
Lx = 2*pi;
Ly = 2*pi;
Nx = 20; % 40; % 6
Ny = 20; % 40; % 5
Nk = 2*(Nx-1)*(Ny-1); % number of triangular elements
dx = Lx/(Nx-1);
dy = Ly/(Ny-1);
nno = 0;
g3 = 1/3; % sqrt(1/3);
rcoef = nquadratic;
%
% Mesh for linear and quadratic functions; nodes and coordinates on corners
%
for jj=1:Ny
  for ii=1:Nx
     nno = nno+1;
     xxn(nno) = (ii-1)*dx;
     yyn(nno) = (jj-1)*dy;
  end
end
%
%
%
if (nquadratic==1)
  %
  % Mesh for quadratic functions; nodes and coordinates on corners
  % 
  % On main rectangle node lines
  %
  for jj=1:Ny
   for ii=1:Nx-1
     nno = nno+1;
     xxn(nno) = 0.5*dx+(ii-1)*dx;
     yyn(nno) = (jj-1)*dy;
   end
  end
  %
  % Between rectangle node lines
  %
  for jj=1:Ny-1
   for ii=1:2*Nx-1
     nno = nno+1;
     xxn(nno) = (ii-1)*0.5*dx;
     yyn(nno) = 0.5*dy+(jj-1)*dy;
   end
  end
  %
end
%
% Index from element and local node number to global node number
%
if (nquadratic == 0)
  Nloc = 3;
  % linear basis functions on triangle 3 nodes
  if (nno ~= Nx*Ny)
     fprintf('Warning: number of nodes nno not equal Nx*Ny in case linear basis functions.\n');
  end
  Index = zeros(Nloc,Nk);
else
  Nloc = 6;
  % quadratic basis functions on triangle 6 nodes
  if ( nno ~= (Nx*Ny+Ny*(Nx-1)+(Ny-1)*(2*Nx-1)) )
     fprintf('Warning: number of nodes nno not equal Nx*Ny,\n in case quadratic basis functions.\n');
  end
  Index = zeros(Nloc,Nk);
end
%
% Number of nodes
%
Nn = nno;
%
% Simple boundary treatment: if node at boundary then last number in Index negative otherwise positive
%
tel = zeros(Nk,1);
for nk=1:Nk
  if (fix(nk/2)*2-nk < 0)
    mk = nk+1; % nk odd
  else
    mk = nk; % nk even
  end
  %
  % Find coordinates of left-bottom corner of the element
  %
  jj = fix((nk-1)/(2*(Nx-1)))+1;
  ii = mk/2-(jj-1)*(Nx-1) ; % Note nominator is always even
  tel(nk,1) = nk; tel(nk,2) = ii; tel(nk,3) = jj;
  if (fix(nk/2)*2-nk < 0) % nk odd
     % local node 1 to 3 triangle oriented counterclockwise; lower right triangle on quad
     Index(1,nk) = (jj-1)*Nx+ii; % global node number
     Index(2,nk) = (jj-1)*Nx+ii+1; % global node number
     Index(3,nk) = jj*Nx+ii+1; % global node number
     if (nquadratic == 1) % define also node 4 to 6
       Index(4,nk) = Nx*Ny+(jj-1)*(Nx-1)+ii; % on line of quads
       Index(5,nk) = Nx*Ny+Ny*(Nx-1)+(jj-1)*(2*Nx-1)+2*ii+1; % in y-middle on right quad line
       Index(6,nk) = Nx*Ny+Ny*(Nx-1)+(jj-1)*(2*Nx-1)+2*ii;   % in y-middle on diagonal
     end
  else % nk even
     % local node 1 to 3 triangle oriented counterclockwise upper left triangle on quad
     Index(1,nk) = (jj-1)*Nx+ii; % global node number same as other triangle
     Index(2,nk) = jj*Nx+ii+1; % global node number; one y-line higher
     Index(3,nk) = jj*Nx+ii; % global node number
     if (nquadratic == 1) % define also node 4 to 6
       Index(4,nk) = Nx*Ny+Ny*(Nx-1)+(jj-1)*(2*Nx-1)+2*ii; % on middle
       Index(5,nk) = Nx*Ny+jj*(Nx-1)+ii; % on top quad line
       Index(6,nk) = Nx*Ny+Ny*(Nx-1)+(jj-1)*(2*Nx-1)+2*ii-1; % on middle
     end
  end
end
%
% Location of boundaries 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
%
switch nquadratic
case 0
   ntfa = (Nx-2)*(Ny-2); % interior nodes
   mnindex = zeros(1,ntfa);
   nb = ntfa;
   ni = 0;
   for m=1:Nn
     if ( (1<=m)&(m<=Nx) | ((Ny-1)*Nx+1<=m)&(m<=Nx*Ny) | mod(m-1,Nx)==0&(m<=Nx*(Ny-1)) | mod(m,Nx)==0&(m<Nx*Ny) )
      % nno is a Dirichlet boundary node
        nb = nb+1;
        mnindex(m) = nb;
        xxnm(nb) = xxn(m);
        yynm(nb) = yyn(m);
     else
      % nno is interior node
        ni = ni+1;
        mnindex(m) = ni;
        xxnm(ni) = xxn(m);
        yynm(ni) = yyn(m);
     end
   end
case 1
   ntfa = (2*Nx-3)*(Ny-1)+(2*Nx-3)*(Ny-2); % interior nodes
   mnindex = zeros(1,ntfa);
   nb = ntfa;
   ni = 0;
   for m=1:Nn
     if (  ((1<=m)&(m<=Nx)) | (((Ny-1)*Nx+1<=m)&(m<=Nx*Ny)) | (mod(m-1,Nx)==0)&(m<=Nx*(Ny-1)) | (mod(m,Nx)==0)&(m<Nx*Ny)  | ...
	 ((Nx*Ny+1<=m)&(m<=Nx*Ny+Nx-1))  | ...
	 ((Nx*Ny+(Ny-1)*(Nx-1)+1<=m)&(m<=Nx*Ny+Ny*(Nx-1)))  | ...
	 ((m-(Nx*Ny+Ny*(Nx-1)+1)>=0)&(mod(m-(Nx*Ny+Ny*(Nx-1)+1),2*Nx-1)==0)) | ...
	 ((m-(Nx*Ny+Ny*(Nx-1)+2*Nx-1)>=0)&(mod(m-(Nx*Ny+Ny*(Nx-1)+2*Nx-1),2*Nx-1)==0)) )
        % m is a Dirichlet boundary node
        nb = nb+1;
        mnindex(m) = nb;
        xxnm(nb) = xxn(m);
        yynm(nb) = yyn(m);
     else
        % m is interior node
        ni = ni+1;
        mnindex(m) = ni;
        xxnm(ni) = xxn(m);
        yynm(ni) = yyn(m);
     end
   end
end
%
% Plot nodes
%
figure(111); clf;
plot(xxn,yyn,'xr','linewidth', 3); hold on;
Lb = -0.1;
axis([Lb 2*Lx Lb 2*Ly]);
xlabel('x');
ylabel('y');