%PDE.m  6.4.04
%%%%%%%%%%%%%%%%%

%Harj. 6
%%%%%%%%%

%Teht. 4
%%%%%%%%%%


>> help delsq

 DELSQ  Construct five-point finite difference Laplacian.
    delsq(G) is the sparse form of the two-dimensional,
    5-point discrete negative Laplacian on the grid G.
    The grid G can be generated by NUMGRID or NESTED.
 
    See also NUMGRID, DEL2, DELSQDEMO.

    
    NUMGRID Number the grid points in a two dimensional region.
    G = NUMGRID('R',n) numbers the points on an n-by-n grid in
    the subregion of -1<=x<=1 and -1<=y<=1 determined by 'R'.
    SPY(NUMGRID('R',n)) plots the points.
    DELSQ(NUMGRID('R',n)) generates the 5-point discrete Laplacian.
    The regions currently available are:
       'S' - the entire square.
       'L' - the L-shaped domain made from 3/4 of the entire square.
       'C' - like the 'L', but with a quarter circle in the 4-th square.
       'D' - the unit disc.
       'A' - an annulus.
       'H' - a heart-shaped cardioid.
       'B' - the exterior of a "Butterfly".
       'N' - a nested dissection ordering of the square.
    To add other regions, edit toolbox/matlab/demos/numgrid.m.
 
    >> G=numgrid('S',5)
G =
     0     0     0     0     0
     0     1     4     7     0
     0     2     5     8     0
     0     3     6     9     0
     0     0     0     0     0

>> full(delsq(G))

ans =

     4    -1     0    -1     0     0     0     0     0
    -1     4    -1     0    -1     0     0     0     0
     0    -1     4     0     0    -1     0     0     0
    -1     0     0     4    -1     0    -1     0     0
     0    -1     0    -1     4    -1     0    -1     0
     0     0    -1     0    -1     4     0     0    -1
     0     0     0    -1     0     0     4    -1     0
     0     0     0     0    -1     0    -1     4    -1
     0     0     0     0     0    -1     0    -1     4

     
     m=3;
     L=numgrid('L',2*m+1)
     L =
     0     0     0     0     0     0     0
     0     1     3     5     7    12     0
     0     2     4     6     8    13     0
     0     0     0     0     9    14     0
     0     0     0     0    10    15     0
     0     0     0     0    11    16     0
     0     0     0     0     0     0     0
>> 

% Moler s. 6
%%%%%%%%%%%%%%

clear ; close all
xv=[0 0 1 1 -1 -1 0]; yv=[0 -1 -1 1 1 0 0]
fill(xv,yv,'y')

n=4; h=1/(n+1);
[x,y]=meshgrid(-1:h:1);
%help inregion

[in,on] = inregion(x,y,xv,yv);

in

hold on

p = find(in-on);   % Sisäsolmujen indeksit (jonouttettuina)
q = find(on);      % Reunasolmujen indeksit
plot(xv,yv,'-',x(p),y(p),'ko',x(q),y(q),'ro')
%axis([-5 5 -4 4])
p     

% Find-kokeilua:
A=hilb(3)
find(A==A(2,1))       % Palauttaa A(:):n indeksit
[i,j]=find(A==A(2,1)) % Palauttaa A:n indeksivektorit i ja j.

q(1:n)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Numgrid, Delsq
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function D = delsq(G)
%DELSQ  Construct five-point finite difference Laplacian.
%   delsq(G) is the sparse form of the two-dimensional,
%   5-point discrete negative Laplacian on the grid G.
%   The grid G can be generated by NUMGRID or NESTED.
%
%   See also NUMGRID, DEL2, DELSQDEMO.

%   C. Moler, 7-16-91.
%   Copyright 1984-2001 The MathWorks, Inc. 
%   $Revision: 5.6 $  $Date: 2001/04/15 12:03:13 $
% G=numgrid('S',5)
[m,n] = size(G);

% Indices of interior points
p = find(G);

% Connect interior points to themselves with 4's.
i = G(p);
j = G(p);
s = 4*ones(size(p));

% for k = north, east, south, west
for k = [-1 m 1 -m]
   % Possible neighbors in k-th direction
%k=-1;
   Q = G(p+k);
   % Index of points with interior neighbors
   q = find(Q);
   % Connect interior points to neighbors with -1's.
   i = [i; G(p(q))];
   j = [j; Q(q)];
   s = [s; -ones(length(q),1)];
end
D = sparse(i,j,s);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Käydään läpi delsq vaiheittain:

clear
G=numgrid('S',5)
[m,n] = size(G);
h=1/(n-1);
% Indices of interior points
p = find(G)
% Connect interior points to themselves with 4's.
i = G(p);j = G(p)
s = 4*ones(size(p));
n2=n-2;

% Oikean puolen b-vektorin kokoaminen:
   % Annettu matriisi F: F(i,j)=f(x_i,y_j).
% Alustus:
F=-0*ones(n2,n2)
ib=i;
b=h^2*F(ib);




% for k = north, east, south, west
for k = [-1 m 1 -m]
   % Possible neighbors in k-th direction
%k=-1;
%k=m;
% k=1;
% k=-m;
   Q = G(p+k);
 %  [i Q]
   % Index of points with interior neighbors
   q = find(Q);
   % Connect interior points to neighbors with -1's.
%   i', G(p(q))'
   i = [i; G(p(q))];
   j', Q(q)'
   j = [j; Q(q)];
%   s', -ones(length(q),1)'
   s = [s; -ones(length(q),1)];
% Oikean puolen b-vektorin kokoaminen:
% Alustettiin ennen for-loopia: b=h^2*F(i);
   qr=find...
   b(qr)=b(qr)+...;
   
end
D = sparse(i,j,s);

%%%%%%%%%%%%%%%%%%%%%%%%%%%

