Finding Grid Point Neighbors Matt Craven for Dr

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(a 7); a 10 = [i, j];

Storing the data in our structures Now that we have a 10, we can use a simple algorithm to find the neighbors of each point, and store them in our structures.

nbfind. m [i, j] = size(a 10); nb = zeros(i, 8); nbor = zeros(i, 8); directions = [ 0, -delta. Y; %North (1) -delta. X, -delta. Y; %Northwest (2) -delta. X, 0; %West(3) -delta. X, delta. Y; %Southwest (4) 0, delta. Y; %South (5) delta. X, delta. Y; delta. X, 0; delta. X, -delta. Y; ]; %Southeast (6) %East (7) %Northeast (8)

nbfind. m, continued for k = 1: i, x = a 10(k, 2); y = a 10(k, 1); number = 0; for dir = 1: 8. nx = x + directions(dir, 1); ny = y + directions(dir, 2); if(ny > 0 & nx > 0) if(a 7(ny, nx)) for c = 1: i, if(a 10(c, 2) == nx) if(a 10(c, 1) == ny) n = c; break; end end number = number + 1; nbnr(k, 1) = number; nb(k, number) = n; nbor(k, number) = dir; 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 a 10. 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 nbor(10, : ) (neighbor order list) 1 6 7 8 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: 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: Summing these, we get:

poisson. m w = zeros(1532, 1532); rhs = zeros(1532, 1); delta. Xsq = 1/delta. X^2; delta. Ysq = 1/delta. Y^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) = delta. Xsq; w(k, S) = delta. Xsq; w(k, E) = delta. Ysq; w(k, W) = delta. Ysq; w(k, k) = -2*(delta. Ysq + delta. Xsq); rhs(k, 1) = 1; else w(k, k) = 1; rhs(k, 1) = 5; end m=max(a 10(: , 1)); n=max(a 10(: , 2)); u = wrhs; w 1=sparse(a 10(: , 1)/delta. Y, a 10(: , 2)/delta. X, u); indxx=delta. X: n; indxy=delta. Y: m; surf(indxy, indxx, w 1');