Transcript Finding Grid Point Neighbors

Finding Grid Point Neighbors
Matt Craven
for:
Dr. Anitescu
Math 1110
Spring 2002
Finding Grid Neighbors
"
"
In order to solve some discretized partial
differerential equations, we need to know the
immediate neighbors of all the points in the
domain.
First, we must choose an appropriate data
structure.
Choosing a data structure
"
In our data structure, we need to maintain three
pieces of information:
–
a vector containing the number of neighbors a node
has
–
a matrix containing an adjacency list of a node's
neighbors
–
another matrix for which entry in the adjacency list is
which neighbor of that node
How the data is stored in the
structures
"
Since we only need to consider the points within
our domain, we create a matrix containing the
coordinates of our points within the domain.
[i,j] = find(a7);
a10 = [i,j];
Storing the data in our structures
Now that we have a10, we can use a simple
algorithm to find the neighbors of each point, and
store them in our structures.
nbfind.m
[i,j] = size(a10);
nb = zeros(i,8);
nbor = zeros(i,8);
directions = [ 0,
%South
-deltaY;
%North (1)
-deltaX, -deltaY;
%Northwest (2)
-deltaX, 0;
%West (3)
-deltaX, deltaY;
%Southwest (4)
0,
deltaY;
(5)
deltaX, deltaY;
deltaX, 0;
deltaX, -deltaY;];
%Southeast (6)
%East (7)
%Northeast (8)
nbfind.m, continued
for k = 1:i,
x = a10(k,2);
y = a10(k,1);
number = 0;
for dir = 1:8.
nx = x + directions(dir, 1);
ny = y + directions(dir, 2);
if(ny > 0 & nx > 0)
if(a7(ny,nx))
for c = 1:i,
if(a10(c, 2) == nx)
if(a10(c,1) == ny)
n = c;
break;
end
end
end
number = number + 1;
nbnr(k,1) = number;
nb(k,number) = n;
nbor(k,number) = dir;
end
end
end
end
How nbfind works
"
"
"
nbfind looks at each point in the domain, and checks each
of the 8 immediate directions to see if a neighbor is
present.
If a neighboring point is in the domain, nbfind looks for
its name in a10.
When the point's name is found, that name is inserted into
the adjacency list, the direction of that point relative to
the current point is placed in the order list, and the
number of neighbors is increased by one.
How the data is stored in each
structure
"
"
"
The number of neighbors (nrnb) vector simply
holds how many points are adjacent to each point
in the domain.
The neighbor adjacency list (nb) holds the name
of the each neighbor, in the order it was found
The neighbor order (nbor) matrix holds a code
describing in which direction the corresponding
neighbor lies, with respect to the node in the
center.
Example entry in our data structure
nbnr(10,:) (number of neighbors)
4
nb(10,:) (neighbor adjacency list)
9
31 30 29 0
0
0
0
nbor(10,:) (neighbor order list)
1
6
7
8
0
0
0
0
9
29
10
30
31
Using this structure to solve the
Poisson equation
"
"
"
Now that we have the data from our domain into
these structures, we can use this information to
solve some discretized differential equations.
To solve these equations, we need to know the
immediate neighbors of every point.
As an example, we will now use this data to
model a solution to the Poisson electrical
potential equation
The Poisson Equation
"
The Poisson equation models electrical potential,
and is given by the partial differential equation:
2
2
 u  u
2 
2 1
x
y
with the boundary condition u = 5, where u(li) is the
potential at point li.
Discretized Poisson Equation
At point C (with neighbors in the four cardinal directions, the
discretized equation has two components:
 u uW  uE  2uC  u uN  uS  2uC


2
2
2
2
x C
(deltaX ) y C
(deltaY )
2
2
Summing these, we get:
uN
uS
uE
uW
2 
2 
2 
2  4uC  1
y
y
x
x
poisson.m
w = zeros(1532, 1532);
rhs = zeros(1532,1);
deltaXsq = 1/deltaX^2;
deltaYsq = 1/deltaY^2;
for k = 1:1532
N = nb(k,find((nbor(k,:) == 1)));
W = nb(k,find((nbor(k,:) == 3)));
S = nb(k,find((nbor(k,:) == 5)));
E = nb(k,find((nbor(k,:) == 7)));
if(~isempty(N) & ~isempty(W) & ~isempty(S) & ~isempty(E))
w(k,N) = deltaXsq;
w(k,S) = deltaXsq;
w(k,E) = deltaYsq;
w(k,W) = deltaYsq;
w(k,k) = -2*(deltaYsq + deltaXsq);
rhs(k,1) = 1;
else
w(k,k) = 1;
rhs(k,1) = 5;
end
end
m=max(a10(:,1));
n=max(a10(:,2));
u = w\rhs;
w1=sparse(a10(:,1)/deltaY,a10(:,2)/deltaX,u);
indxx=deltaX:deltaX:n;
indxy=deltaY:deltaY:m;
surf(indxy,indxx,w1');