clc
clear
%%--------------User Inputs - Pre processing-------------------------------
order = 3; % Gauss quadrature order
pol_degree = 2; % Lagrange polynomial degree of shape functions
nelem1 = 5; %%Number of elements for x coord
nelem2 = 5; %%Number of elements for y coord
nelem = nelem1 * nelem2; % Total Number of elements
nelem1;

NN1 = nelem1 * pol_degree + 1; %% Number of Nodes x coord
NN2 = nelem2 * pol_degree + 1; %% Number of Nodes y coord
NN = NN1 * NN2; %% Total # Nodes
domain1 = 1.0 ;% The actual domain of the ODE 0<x<1
domain2 = 1.0 ;

%%Provide node coordinates given the number of elements, pol_degree and domain
%%*******MESH**********************
% Equally spaced quadratic (pol_degree=2) elements
%%domain1 = x
%%domain2 = y
dx = domain1/nelem1/pol_degree; % x distance between x nodes
dy = domain2/nelem2/pol_degree; % y distance between y nodes

xnode(1,1) = 0;
xnode(2,1) = 0;
%% Define X-node Coordinates

for j = 1 : NN2
for i = 1 : NN1
xnode(1,i + (j-1) * NN1) = xnode(1,1) + (i-1) * dx;
xnode(2,i + (j-1) * NN1) = xnode(2,1) + (j-1) * dy;

end
end
% figure
% scatter((xnode(1,:)), xnode(2,:));
% title('Mesh Plot')
% hold on
% axis equal
% grid on

%% Way we measure nodes in 2D is counter clockwise on the exterior nodes
%% then interior nodes.
%%Provide gcon == connectivity matrix

gcon = zeros((pol_degree+1)^2, nelem);
for j = 1:nelem2
for i = 1:nelem1
k = nelem1 * (j-1) + i;

%% 2 * nelem + 1 = how we increase from row 1 nodes to row 2
gcon(1, k) = 1 + pol_degree*(i-1) + (j-1)*pol_degree*NN1;
gcon(2, k) = 3 + pol_degree*(i-1) + (j-1)*pol_degree*NN1;
gcon(3, k) = 3 + pol_degree*(i-1) + (j-1)*pol_degree*NN1 +...
(2*nelem1 + 1)*2;

gcon(4, k) = 1 + pol_degree*(i-1) + (j-1)*pol_degree*NN1 +...
(2*nelem1 + 1)*2;
gcon(5, k) = 2 + pol_degree*(i-1) + (j-1)*pol_degree*NN1;
gcon(6, k) = 1 + 3 + pol_degree*(i-1) + (j-1)*pol_degree*NN1 +...
2*nelem1;
gcon(7, k) = 2 + pol_degree*(i-1) + (j-1)*pol_degree*NN1 +...
(2*nelem1 + 1)*2;
gcon(8, k) = 2 +pol_degree*(i-1) + (j-1)*pol_degree*NN1 +...
2*nelem1 ;
gcon(9, k) = 3 + pol_degree*(i-1) + (j-1)*pol_degree*NN1 +...
2*nelem1;
end

Introduction to Finite Element Coding Homework 2D Consider the Po...