% Compute the L shaped region problem using quadralateral elements with
% no quadrature.  The L shaped region is defined on a 2 unit by 2 unit
% square, with coordinates running from -1 to 1 on each axis.  The lower
% left quadrant (x1<0, x2<0) is blanked out, leaving a L shaped
% comutational domain.  k*Laplacian(q) - s = 0 is solved, with k=1 and
% s=1 on the interior.  A Dirichlet condition (q=0) is imposed along the
% entire boundary.

% Declare global variables.  These variables must be made available as
% global variables for the assembly routine to work.

elwide = 10;				% number of elements along the
                                        % width of the grid
eltall = 10;				% number of elements along the
                                        % height of the grid
					
% Construct the grid seeds.  Be sure that X1seed and X2seed both include
% 0.  
X1seed=linspace(-1,1,elwide+1)
X2seed=linspace(1,-1,eltall+1)

% construct the nodal coordinate lists
[X1 X2] = meshgrid(X1seed, X2seed);
nnodes = prod(size(X1));

% construct the element connectivity table
elcon = quadcon(elwide,eltall);

% Scan through all nodes and select the ones for which Dirichlet
% conditions will be enforced. Also, select the nodes outside the
% computational domain for blanking.  This will be enforced in exactly
% the same manner as Dirichlet B.C's.
diri_list = do_Lshape2D_dlist(X1,X2);

% calculate d(etai)/d(xj) and determinants for each element
[DETe,DETeI,ETAe]=do_ETAe_q0(X1,X2,elcon);

% calculate diffusion matrix for each element
DIFFe = do_DIFFe_q0(ETAe);

% assemble the system diffusion matrix
DIFF = asjac([],[],[],-1,DIFFe,[],...
	     elcon,DETe,nnodes);

% assemble the residual
B200L = 1/9 * [4 2 1 2;...
	       2 4 2 1;...
	       1 2 4 2;...
	       2 1 2 4];
SRC = ones(size(X1));
b = asres([],[],[],1,B200L,SRC,...
	  elcon,DETe,nnodes);

% apply the BC's/blanking
[DIFF_diri b_diri] = do_diri(DIFF,b,diri_list);

% solve the linear system
Q = DIFF_diri\b_diri;

% reshape Q from a column matrix to an array of the same shape as X1 for
% plotting purposes.
Q_filled = zeros(size(X1));
Q_filled(:) = Q;

surfc(X1,X2,Q_filled);
