# sandbox/easystab/dif2D.m

In this code, we build the 2D differentiation matrices of a rectangular domain.

To learn how to use this function, please see diffmat_2D_dif2D.m.

function [D,l,X,Y,Z,I]=dif2D(d,x,y)

Nx=length(x);
Ny=length(y);
NN=Nx*Ny;


# 2D differentiation matrices

To build the 2D differentiation matrices from the 1D ones, we use the kron operator. Please see diffmat_2D.m#d-matrices-1.

% differentiation
D.x=kron(d.x,speye(Ny));
D.xx=kron(d.xx,speye(Ny));
D.y=kron(speye(Nx),d.y);
D.yy=kron(speye(Nx),d.yy);


# Integration weights in 2D

% integration
D.w=d.wy(:)*d.wx(:)';
D.wx=ones(Ny,1)*d.wx(:)';
D.wy=d.wy(:)*ones(1,Nx);

D.w=reshape(D.w,Ny,Nx);
D.wx=reshape(D.wx,Ny,Nx);
D.wy=reshape(D.wy,Ny,Nx);


# Location vectors

To impose the boundary conditions we need to know where in the vector representation of the 2D variable are the different boundaries and corners. For this we build here a few location vectors stored as the fields of the structure l. To learn about the use of location vectors, please see pedagogy#location-vectors-and-matrices. To learn how the array variable is transformed into a vector, see diffmat_2D.m; there you will as well understand the structure of the mesh in 2D.

Here we discuss the different fields of l and how they are built. You see below the structure of the mesh: The different boundaries of the 2D mesh

We start by building the long vector of all the indices 1:NN, where NN=NxNy* is the total number of degrees of freedom of the 2D mesh. This vector we transform into the shape of the mesh with the function reshape. We store this in the array dom (standing for “domain”). Once this done, we have the link between the positions in the array and the locations in the vector. Since the indices in the array increase with position, we have that:

• the first line of the array dom corresponds to the bottom boundary cells,
• the last line of the array dom corresponds to the top boundary cells,
• the first column of the array dom corresponds to the left boundary cells,
• the last column of the array dom corresponds to the right boundary cells.
• The corners of the array dom correspond to the corner cells.

This means that when you display the array as Octave/Matlab does on screen, you see your mesh upside-down. To avoid this we could have chosen to have the y position of cells decreasing with index instead of having them increase with index, but I prefered to keep the same convention for x and y.

In this generic function, I chose to store the side boundaries without the corners. Depending on your physical system, you may chose to include the corners in the location vectors for instance in the right and left boundaries, and keep the top and bottom boundaries without corners, see for instance ???

% locations
dom=reshape(1:NN,Ny,Nx);
right=dom(2:Ny-1,end);    l.right=right(:);
left=dom(2:Ny-1,1);       l.left=left(:);
top=dom(end,2:Nx-1);      l.top=top(:);
bot=dom(1,2:Nx-1);    l.bot=bot(:);
cor=dom([1 end],[1 end]);     l.cor=cor(:);

l.cbl=cor(1); % bottom left corner
l.ctl=cor(2); % top left corner
l.cbr=cor(3); % bottom right corner
l.ctr=cor(4); % top right corner


[X,Y]=meshgrid(x,y);

Z=spalloc(NN,NN,0); I=speye(NN);