close all
gf_workspace('clear all');
show_figures = 1;
%
% creiamo una mesh sul quadrato [0,1]x[0,1] di passo hx=hy=h :
h = 0.1;
m = gf_mesh('triangles_grid',[0:h:1],[0:h:1]);
% e facciamone il grafico:
if 0
  figure(1); gf_plot_mesh(m);
  pause
end
%
% applichiamo una formula di integrazione esatta per gli elementi geometrici di "m" :
mim = gf_mesh_im(m, gf_integ('IM_EXACT_SIMPLEX(2)'));
%
% creiamo la struttura dati "mfu" per la soluzione discreta di un campo scalare (dim="1") definito sulla mesh "m" :
mfu = gf_mesh_fem(m,1);
% ed assegnamo elementi finiti di tipo P1 a tutti gli elementi di "mfu" :
gf_mesh_fem_set(mfu,'fem',gf_fem('FEM_PK(2,1)'));
%
% creiamo la struttura dati "mfd" per il termine noto definito sulla mesh "m" :
mfd = gf_mesh_fem(m,1);
% ed assegnamo elementi finiti di tipo P1 a tutti gli elementi di "mfd" :
gf_mesh_fem_set(mfd,'fem',gf_fem('FEM_PK(2,1)'));
N = gf_mesh_fem_get(mfd, 'nbdof');
%
% costruiamo una descrizione del bordo
border = gf_mesh_get(m,'outer faces');
% e segnamolo come segmento n. 1 :
gf_mesh_set(m, 'boundary', 1, border);
if 0
  gf_plot_mesh(m, 'regions', [1]); % i segmenti del bordo sono in rosso.
  pause;
end
%
% definiamo la condizione iniziale "U0" e ne calcoliamo l'interpolante sulla mesh "mf" :
Ud = gf_mesh_fem_get(mfu, 'eval', { '0' });
pid=gf_mesh_fem_get(mfu,'dof from cvid',border);
%
% soluzione senza i "model bricks" (e quindi a piu' basso livello):
%
% condizioni al contorno di Dirichlet:
[H,R] = gf_asm('dirichlet', 1, mim, mfu, mfd, gf_mesh_fem_get(mfd,'eval',1),Ud);
[null,Ud]=gf_spmat_get(H,'dirichlet nullspace', R);
Ud(find(sum(null,2)==0)') = 0.9;
%
% discretizzazione della PDE:  
beta = 0;
M = gf_asm('mass matrix',mim, mfu);
K = gf_asm('laplacian',mim, mfu,mfd,ones(1,N));
%A = K + beta * M;
% elimino le condizioni al contorno di Dirichlet:
F=null'*(-(K+beta*M)*Ud(:));
K=null'*K*null;
M=null'*M*null;
%
Tf = 0.01;
dt = 0.0005;
%
% Analisi della matrice di iterazione:
if 1
  rho_eig = sort(eig(full(inv(M)*K)))
  sort(dt*(beta+rho_eig))
  autov_mat_iteraz = sort(eig((1-dt*beta)*eye(size(M))-dt*inv(M)*K))
end
%
% definiamo la condizione iniziale "U0":
U0 = Ud;
if show_figures
  figure(1); gf_plot(mfu,U0,'zplot','on'); colorbar; title('initial value'); view(2);           
end
%
u_nd = null'*U0';
for i=1:floor(Tf/dt)
  u_nd = M \ (((1-dt*beta)*M-dt*K)*u_nd + dt*F);
  U=null*u_nd+Ud(:);
  if show_figures
    Ufig=gf_compute(mfu,U','interpolate on',mfd);
    figure(1); gf_plot(mfu,Ufig,'zplot','on'); colorbar; title('solution'); view(2); 
    drawnow, pause(0.4);
  end
  disp(['t = ' num2str(i*dt) '   min(U) = ' num2str(min(U))])
end