function demo_poisson
% Interactive demo: solving Poissons equation with the FEMB program collection 

% Initialize data
global FEMB_Data;
FEMB_Data.H=9;
FEMB_Data.n=3;
FEMB_Data.step=0;
FEMB_Data.cellind=[4,4];
FEMB_Data.funstring='@(x,y) sin(2*pi*y)';
FEMB_Data.f=@(x,y) sin(2*pi*y);
clf reset

% Initialize GUI
set(gcf,'units','normalized');
set(gcf,'WindowButtonDownFcn',@FEMB_Button);
set(gcf,'KeyPressFcn',@FEMB_Key);

% Text for title
FEMB_Data.uicontrols(1)=uicontrol('Units','normalized', ...
  	'BackgroundColor',get(gcf,'Color'), ...
	'HorizontalAlignment','center', ...
        'FontWeight','bold',...
        'FontSize',20,...
	'Position',[.0,.90,.99,.09],...
	'String','Poisson Solver with Weighted B-Splines',...
	'Style','text');

% Text for degree
FEMB_Data.uicontrols(2)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.05,.85,.1,.05],...
	'String','Degree:',...
	'Style','text');

% Input for degree
FEMB_Data.uicontrols(3)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','right', ...
	'Position',[.15,.85,.03,.05],...
	'String',sprintf('%1d',FEMB_Data.n),...
	'Style','edit',...
	'Callback',@degree_update);

% Text for grid cells
FEMB_Data.uicontrols(4)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','right', ...
	'Position',[.18,.85,.15,.05],...
	'String','Grid cells:',...
	'Style','text');

% Input for grid cells
FEMB_Data.uicontrols(5)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','right', ...
	'Position',[.35,.85,.05,.05],...
	'String',sprintf('%3d',FEMB_Data.H),...
	'Style','edit',...
	'Callback',@cellnum_update);

% Start/next button
FEMB_Data.uicontrols(6)=uicontrol('Units','normalized', ...
	'BackgroundColor',get(gcf,'Color'), ...
	'String','Start', ...
	'Position',[.85,.85,.1,.05], ...
	'Callback',@nextstep,...
	'Style','pushbutton');

% Text for Function
FEMB_Data.uicontrols(8)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','right', ...
	'Position',[.42,.85,.09,.05],...
	'String','Force:',...
	'Style','text');

% Input for Function
FEMB_Data.uicontrols(9)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.52,.85,.3,.05],...
	'String',sprintf('%s',FEMB_Data.funstring),...
	'Style','edit',...
	'Callback',@function_update);

% Text for user info/subtitle
FEMB_Data.uicontrols(7)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'ForegroundColor',[0,0,0], ...
        'FontWeight','bold',...
        'FontSize',12,...
        'HorizontalAlignment','center', ...
	'Position',[.0,0,.99,.05],...
	'String','',...
	'Style','text');

% Different axes for display
FEMB_Data.axis(1)=subplot('position',[.05,.1,.9,.68]);
FEMB_Data.axis(2)=0;
FEMB_Data.axis(3)=0;
myplot;

%%%%%  Callback Functions 
% mouse button pressed 
function FEMB_Button(varargin)
global FEMB_Data;
if FEMB_Data.step == 6 & gca ==  FEMB_Data.axis(1)
  if strcmp(get(gcf,'SelectionType'),'normal') % Left mouse
    p=get(FEMB_Data.axis(1),'CurrentPoint');
    if p(1) >= 0 & p(1)  <= 1
      if p(3) >= 0 & p(3)  <= 1
	FEMB_Data.cellind = ceil(FEMB_Data.H*p([1,3]));
	myplot;
      end    
    end
  end
end

% key pressed 
function FEMB_Key(varargin)
global FEMB_Data;
c=double(get(gcf,'currentcharacter'));
if FEMB_Data.step == 5 
  upbound=FEMB_Data.H+FEMB_Data.n;
end
if FEMB_Data.step == 6 
  upbound=FEMB_Data.H;
end
if FEMB_Data.step == 5 | FEMB_Data.step == 6 
  switch c
   case 29
    if FEMB_Data.cellind(1) < upbound
      FEMB_Data.cellind(1)=FEMB_Data.cellind(1)+1;
      myplot
    end;
   case 28
    if FEMB_Data.cellind(1) > 1
      FEMB_Data.cellind(1)=FEMB_Data.cellind(1)-1;
      myplot
    end;
   case 30
    if FEMB_Data.cellind(2) < upbound
      FEMB_Data.cellind(2)=FEMB_Data.cellind(2)+1;
      myplot
    end;
   case 31
    if FEMB_Data.cellind(2) > 1
      FEMB_Data.cellind(2)=FEMB_Data.cellind(2)-1;
      myplot
    end;
  end
end

% force input
function function_update(varargin)
global FEMB_Data;
FEMB_Data.funstring = get(FEMB_Data.uicontrols(9),'String');
FEMB_Data.f= eval(FEMB_Data.funstring);
restart;

% degree input
function degree_update(varargin)
global FEMB_Data;
n= str2double(get(FEMB_Data.uicontrols(3),'String'));
if n < 1 | n > 9
  set(FEMB_Data.uicontrols(3),'String',sprintf('%1d',FEMB_Data.n));
else
  FEMB_Data.n =n;
  restart;
end

% grid cells input
function cellnum_update(varargin)
global FEMB_Data;
FEMB_Data.H = str2double(get(FEMB_Data.uicontrols(5),'String'));
restart;


% restart subroutine
function restart(varargin)
global FEMB_Data;
FEMB_Data.step = 0;
set(FEMB_Data.uicontrols(6),'String','Start');
set(FEMB_Data.uicontrols(7),'String','');   
set(FEMB_Data.uicontrols(7),'visible','on');   
if FEMB_Data.axis(3)
  axes(FEMB_Data.axis(3))
  cla reset;
  delete(FEMB_Data.axis(3));
  FEMB_Data.axis(3)=0;
end
FEMB_Data.axis(1)=subplot('position',[.05,.1,.9,.68]);
myplot;

% next button pressed 
function nextstep(varargin)
global FEMB_Data;
if FEMB_Data.step < 10
  FEMB_Data.step =FEMB_Data.step+1;
else
  restart;
end
if FEMB_Data.step > 0  
  set(FEMB_Data.uicontrols(6),'String','Next');
end
if FEMB_Data.step == 4  
  FEMB_Data.PAR=set_par(FEMB_Data.H,FEMB_Data.n);
  FEMB_Data.PAR.info=0;FEMB_Data.PAR.values=15;
  
  [FEMB_Data.int,FEMB_Data.cells]=integrate_2d(@w_shovel,FEMB_Data.H,FEMB_Data.PAR);
end
if FEMB_Data.step == 5
  FEMB_Data.cellind=FEMB_Data.n*[1,1]+1;   
end
if FEMB_Data.step == 6
  FEMB_Data.cellind=[1,1];   
  FEMB_Data.intpts=0;   
end
myplot;

% plot-function
function myplot(varargin)
global FEMB_Data;
  H=FEMB_Data.H;
  n=FEMB_Data.n;
  wf=@w_shovel;
  f=FEMB_Data.f;
  rotate3d off;

  % date for grid and boundary
  [x,y]=ndgrid([0:.01:1]);
  Xg = [[0:H]/H, zeros(1,H+1); [0:H]/H, ones(1,H+1)];
  Yg = [zeros(1,H+1), [0:H]/H; ones(1,H+1), [0:H]/H];
  
switch FEMB_Data.step 
 case 0 % start
  cla reset;
  set(FEMB_Data.uicontrols(7),'String','');
  axis off
 case 1 % domain
  cla reset;
  [x,y]=ndgrid([0:.01:1]);
  contourf(x,y,w_shovel(x,y),[0,inf]); 
  axis equal;axis off;
  colormap winter;
  set(FEMB_Data.uicontrols(7),'String','Domain');  
 case 2 % weight function
  cla reset;
  w=wf(x,y);w(w<0)=0;
  surf(x,y,w);hold on;
  contour(x,y,wf(x,y),[0 0],'-k','linewidth',2);
  shading interp;colormap winter
  view(3);rotate3d on;
  axis on;
%  axis equal;axis off
  set(FEMB_Data.uicontrols(7),'String','Weight function');  
 case 3 % grid
  cla reset;
  contour(x,y,wf(x,y),[0 0],'-k','linewidth',2);hold on
  plot(Xg,Yg,'-k')
  axis([0 1 0 1]);axis equal;axis off;
  set(FEMB_Data.uicontrols(7),'String','Grid');  
 case 4 %  cell classification
  cla reset;
  hold on
  plot(Xg,Yg,'-k')
  axis([0 1 0 1]);axis equal;axis off;
  color={[0,0,1],[0,.5,.75],[0,1,.5]};
  for kx=1:H,for ky=1:H
      fill((kx-[1,0,0,1])/H,(ky-[0,0,1,1])/H,color{FEMB_Data.cells(kx,ky)+2},...
	   	   'facealpha',.4);
  end;end
  contour(x,y,wf(x,y),[0 0],'-k','linewidth',2);
  set(FEMB_Data.uicontrols(7),'String','Cell classification');  
 case 5 % B-splines
  cla reset;
  contour(x,y,wf(x,y),[0 0],'-k','linewidth',2);hold on
  plot(Xg,Yg,'-k')
  U=zeros(H+n);
  U(FEMB_Data.cellind(1),FEMB_Data.cellind(2))=1;
  [uXY,X,Y] = spline_2d(U,H,n,25);
  uXY(uXY==0)=NaN;
  w=wf(X,Y);
  w(w<0)=0;
  zm=max(max(uXY.*w));
  surf(X,Y,uXY.*w,'edgecolor','none');
  hold off;
  view(-30,60)
  axis([0 1 0 1 0 zm+.1]) ;axis off;
  set(FEMB_Data.uicontrols(7),'String',...
    sprintf('B-spline (%d,%d), use arrow keys to move',FEMB_Data.cellind));  
 case 6 % integration points
  if FEMB_Data.intpts == 0
    FEMB_Data.axis(1)=subplot('position',[.05,.1,.4,.68]);cla reset;
    contour(x,y,wf(x,y),[0 0],'-k','linewidth',2);hold on
    plot(Xg,Yg,'-k')
    axis equal ;axis off
    color={[0,0,1],[0,.5,.75],[0,1,.5]};
    for kx=1:H,for ky=1:H
	fill((kx-[1,0,0,1])/H,(ky-[0,0,1,1])/H,color{FEMB_Data.cells(kx,ky)+2},...
	     'facealpha',.4);      
      end;end    
      FEMB_Data.intpointbox=...
	  fill((FEMB_Data.cellind(1)-[1,0,0,1])/H,(FEMB_Data.cellind(2)-[0,0,1,1])/H, .1*[1,1,1]);
      axis([0,1,0,1]);
      FEMB_Data.axis(2)=subplot('position',[.5,.1,.4,.68]);
      contour(x,y,wf(x,y),[0 0],'-k','linewidth',2);hold on
      plot(Xg,Yg,'-k')
      p=FEMB_Data.int.p{FEMB_Data.cellind(1),FEMB_Data.cellind(2)};  
      FEMB_Data.intpointh=plot(p(:,1),p(:,2),'b.','Markersize',15);
      axis equal ;axis off
      axis([FEMB_Data.cellind(1)+[-1.2,.2],FEMB_Data.cellind(2)+[-1.2,.2]]/FEMB_Data.H);
      hold off;
      set(FEMB_Data.uicontrols(7),'String','Gauss points, click on cell or use arrow keys to move');  
      FEMB_Data.intpts =1;
  else
    axes(FEMB_Data.axis(2));
    set(FEMB_Data.intpointbox,'Xdata',(FEMB_Data.cellind(1)-[1,0,0,1])/H);
    set(FEMB_Data.intpointbox,'Ydata',(FEMB_Data.cellind(2)-[0,0,1,1])/H);
    p=FEMB_Data.int.p{FEMB_Data.cellind(1),FEMB_Data.cellind(2)};  
    set(FEMB_Data.intpointh,'Xdata',p(:,1));
    set(FEMB_Data.intpointh,'Ydata',p(:,2));
    axis([FEMB_Data.cellind(1)+[-1.2,.2],FEMB_Data.cellind(2)+[-1.2,.2]]/FEMB_Data.H);
  end
 case 7 % assembling
  cla reset;
  FEMB_Data.axis(1)=subplot('position',[.05,.1,.9,.68]);
  contour(x,y,wf(x,y),[0 0],'-k','linewidth',2);hold on
  plot(Xg,Yg,'-k')
  axis([0 1 0 1]), axis equal ;axis off  
  set(FEMB_Data.uicontrols(6),'visible','off');     
  if H < 21
    set(FEMB_Data.uicontrols(7),'String',...
      sprintf('Assembling Ritz-Galerkin system, speed reduction %5.3fs per cell',1/H/2));
  else
    set(FEMB_Data.uicontrols(7),'String',...
      sprintf('Assembling Ritz-Galerkin system, speed reduction %5.3fs per column',1/H));
  end    
  [FEMB_Data.G,FEMB_Data.F] = myassemble_2d(@p_poisson,@q_poisson,...
					  f,wf,...
					  H,n,FEMB_Data.int,FEMB_Data.PAR);
  set(FEMB_Data.uicontrols(7),'String','Ritz-Galerkin system assembled'); 
  set(FEMB_Data.uicontrols(6),'visible','on');   
 case 8 % solving
  cla reset;
  set(FEMB_Data.uicontrols(7),'String','Solving system with cg iteration');    
  set(FEMB_Data.uicontrols(6),'visible','off');   
  FEMB_Data.U = mysolve_2d(FEMB_Data.G,FEMB_Data.F,FEMB_Data.PAR);  
  set(FEMB_Data.uicontrols(7),'String','Ritz-Galerkin system solved');    
  set(FEMB_Data.uicontrols(6),'visible','on');   
  axis off
 case 9 % solution
  cla reset;
  [FEMB_Data.sol,FEMB_Data.X,FEMB_Data.Y,FEMB_Data.R,FEMB_Data.e] =...
      evaluate_2d(@p_poisson,@q_poisson,f,wf, ...
		  wf,FEMB_Data.U,FEMB_Data.H,n,FEMB_Data.PAR);  
  w=wf(FEMB_Data.X,FEMB_Data.Y);
  FEMB_Data.sol(w<0)=0;
  contour(FEMB_Data.X,FEMB_Data.Y,w,[0 0],'-k','linewidth',2);hold on;
  surf(FEMB_Data.X,FEMB_Data.Y,FEMB_Data.sol);
  shading interp;colormap winter
  view(3);rotate3d on;  axis on
  set(FEMB_Data.uicontrols(7),'String','Solution');    
 case 10 % error
  cla reset;
  if n > 1
    surf(FEMB_Data.X,FEMB_Data.Y,FEMB_Data.R);
    shading interp;colormap hsv
    view(3);rotate3d on;  axis on
  else
    text(.5,.5,'Residual not available for degree 1', ...
	 'fontsize',20,'horizontalalignment','center');
    axis off
  end
  set(FEMB_Data.uicontrols(7),'String','');
  set(FEMB_Data.uicontrols(7),'visible','off');
  FEMB_Data.axis(3)=subplot('position',[.0,0,.99,.05]);
  axis([0,1,0,1]);axis off
  text(.5,.5,sprintf('Relative L_2 error: || \\Delta u -\\Delta u_h || / || \\Delta u || =  %5g',FEMB_Data.e),...
       'FontWeight','bold',...
       'FontSize',12,...
       'interpreter','tex',...
       'HorizontalAlignment','center');
  set(FEMB_Data.uicontrols(6),'String','Restart');  
end



%%%%% Modified FEMB functions (plot during computation)

function [G,F] = myassemble_2d(p,q,f,w,H,n,int,PAR)
% [G,F] = assemble_2d(p,q,f,w,H,n,int,PAR)
% assembly of the Ritz-Galerkin system

% initialization 
G = zeros(H+n,H+n,2*n+1,2*n+1); F = zeros(H+n,H+n);

% loop over grid cells 
work = 0; 
for mx = 1:H
  if H > 20
    fill((mx+[-1,0,0,-1])/H,([0,0,H,H])/H,0.*[1,1,1],'facealpha',.6);
    pause(1/H)
  end    
  for my = 1:H
  work = work + 100/H^2;
  if H< 21
    fill((mx+[-1,0,0,-1])/H,(my+[-1,-1,0,0])/H,0.*[1,1,1],'facealpha',.6);
    pause(1/H/2)
  end
  area = sum(int.c{mx,my}); 
  if area > 0
    % values of relevant weighted b-splines for cell m
    px = int.p{mx,my}(:,1); py = int.p{mx,my}(:,2); 
    [B,Bx,By] = bspline_2d(w,n,H,px,py); 
    Pxy = p(px,py); Qxy = q(px,py); Fxy = f(px,py);
 
    % i: indices of b-splines overlapping cell m
    % b-spline i has segment n+1-(i-m) on cell m 
    for ix = mx:mx+n, for iy = my:my+n

      % linear form
      Ix=n+1-ix+mx; Iy=n+1-iy+my;
      F(ix,iy) = F(ix,iy)+int.c{mx,my}'*(Fxy.*B{Ix,Iy});
      for jx=mx:mx+n, for jy=my:my+n

        % bilinear form
        Jx=n+1-jx+mx; Jy=n+1-jy+my;
        G(ix,iy,jx-ix+n+1,jy-iy+n+1) = ... 
          G(ix,iy,jx-ix+n+1,jy-iy+n+1)+int.c{mx,my}'* ... 
          bilinear(Pxy,Qxy,B{Ix,Iy},Bx{Ix,Iy},By{Ix,Iy}, ... 
            B{Jx,Jy},Bx{Jx,Jy},By{Jx,Jy});
      end, end
    end, end
  end
end, end

% diagonal entries for irrelevant b-splines -> 1
% corresponding entries of right side -> 0
maxG = norm(G(:),inf);
diagG = G(:,:,n+1,n+1); 
F(abs(diagG)<=maxG*PAR.tol_zero^2) = 0;
diagG(abs(diagG)<=maxG*PAR.tol_zero^2) = 1; 
G(:,:,n+1,n+1) = diagG;

function U = mysolve_2d(G,F,PAR)
% U = solve_2d(G,F,PAR)
% conjugate gradient solution of the Ritz-Galerkin system

% initialization
n = (size(G,3)-1)/2;
H = size(G,1)-n; 
norm_F = norm(F(:));

% symmetric diagonal preconditioning
DG = ones(H+3*n,H+3*n); 
ind = [1:H+n];
DG(n+ind,n+ind) = 1./sqrt(squeeze(G(:,:,n+1,n+1)));
for sx=1:2*n+1, for sy=1:2*n+1
  G(:,:,sx,sy) = G(:,:,sx,sy).*DG(n+ind,n+ind).* ... 
    DG(sx-1+ind,sy-1+ind);
end, end
F = F.*DG(n+ind,n+ind);

% start values
U = zeros(H+n,H+n);
D = -F; V = -D; W = 0*D;
gamma = D(:)'*D(:); 
tol = gamma*PAR.tol_solve^2;
hold on;axis equal
fill([tol/10,1,1,tol/10],[0,0,1,1],.75*[1,1,1]);
h1=fill([tol/10,1,1,tol/10],[0,0,1,1],.75*[0,1,0]);
set(gca,'XScale','log');
axis([tol/10,10,-.3,1.3]);axis off
plot([tol,tol],[-.5,1.5],'-k','linewidth',2);
text(tol,1.5,sprintf('tol=%5g',tol));
%title(sprintf('Speed reduction %6.3gs per step',1/H),'fontweight','bold');
h=text(exp((log(tol)+log(10))/2),2.3,...
    sprintf('Speed reduction %6.3gs per step',1/H),...
     'fontweight','bold',...
     'horizontalalignment','center');
h2=text(gamma,-.5,sprintf('error=%5g',gamma));

% iteration 
for steps = 1:PAR.iterations
  if PAR.info & mod(steps,10)==0, fprintf('\b\b\b\b\b\b%6d',steps), end      
   W = mult(G,V,H,n);
   alpha = gamma/(V(:)'*W(:)); beta = 1/gamma;
   U = U + alpha*V;
   D = D + alpha*W; 
   gamma = D(:)'*D(:);
   set(h1,'Xdata',[tol/10,gamma,gamma,tol/10]);pause(1/H);
   set(h2,'String',sprintf('error=%5g',gamma));
   set(h2,'position',[gamma,-.5])
   if gamma < tol
     break
   end
   beta = beta*gamma; 
   V = -D + beta*V;
end 
delete(h);

% error 
error = (mult(G,U,H,n)-F)./DG(n+ind,n+ind);

% reverse preconditioning 
U = U.*DG(n+ind,n+ind);

%%%%% local FEMB functions 

function [B,Bx,By] = bspline_2d(w,n,H,px,py)
% values and derivatives of weighted b-splines

% weight function
[W,Wx,Wy] = w(px,py);

% b-splines
px = mod(px*H,1); py = mod(py*H,1); 
[sx,dsx] = b_spline(px,n); [sy,dsy] = b_spline(py,n);
for kx=1:n+1, for ky=1:n+1
  b{kx,ky} = sx(:,kx).*sy(:,ky);
  bx{kx,ky} = dsx(:,kx).*sy(:,ky)*H;
  by{kx,ky} = sx(:,kx).*dsy(:,ky)*H;
end, end

% weighted b-splines for relevant segments
for kx=1:n+1, for ky=1:n+1
  B{kx,ky} = W.*b{kx,ky};
  Bx{kx,ky} = W.*bx{kx,ky}+Wx.*b{kx,ky};
  By{kx,ky} = W.*by{kx,ky}+Wy.*b{kx,ky};
end, end

function GIJ = bilinear(Pxy,Qxy,BI,BIx,BIy,BJ,BJx,BJy)
% values of the bilinear form

GIJ = Pxy.*(BIx.*BJx+BIy.*BJy) + Qxy.*BI.*BJ;

function W = mult(G,V,H,n)
% sparse matrix/vector multiplication

V0 = zeros(H+3*n,H+3*n); 
V0(n+1:H+2*n,n+1:H+2*n) = V;
W = zeros(H+n,H+n);
for sx=1:2*n+1, for sy=1:2*n+1
  W = W + G(:,:,sx,sy).*V0(sx:end-2*n-1+sx,sy:end-2*n-1+sy); 
end, end







