function [u,p,e,t]=femsol(g,b,c,a,f,hmax)
% [u,p,e,t]=femsol(g,b,c,a,f,hmax)
%
% Yhtalona Poisson: div(c*grad u)+a u = f
%
% Kirjoita parametrifunktiot c,a ja f (R^2 -> R) hipsuissa,
% joko tyyliin: '2' tai 'x.^2 + y.^2' tai 'foo(x,y)'.  
%
% Muodosta tehtavan geometria g ja reunaehdot b  pdetoolin avulla. 
% Sallitut reunaehdot:
%   u = r      Dirichlet  ja
%   u' = g     Neumann.
%
% hmax on maksimi kolmion sivun pituus, kokeile piirtamalla:
% pdeplot(p,e,t) 

% kolmiointi: p,e,t (pdetool)
[p,e,t] = initmesh(g,'Hmax',hmax,'jiggle','off');

% massa- ja jäykkyysmatriisit sekä oikean puolen vektori F
[K,M,F]=lassema(p,t,c,a,f);

% reunaehtojen kontribuutiot, Dirichlet: R,H,  Neumann: G
[G,H,R]=lassemb(b,p,e);

% ratkaiseminen               
[tmp idp] = find(H);       % Dirichlet-pisteiden indeksit
tmp=ones(size(H,2),1); tmp(idp)=0;
[inp] = find(tmp);         % Neumann-pisteiden indeksit

aa = K+M;  aa(idp,:)=H; 
bb = F+G;  bb(idp,:)=R;

u=aa\bb;


%   Author: Reija Ristilä <reija@hut.fi>
%   Last modified: 1999-05-11 12:11:17 EEST

