function [Operators] = Q1FEM(xgl)
% function [Operators] = Q1FEM(x)
%
% Based on grid x, will create a mass and stiffness
% matrix for use in elliptic problems.
%
% Builds tensor product like Q1 FEM operators
% This is mostly used for preconditioning purposes
%
% % Amik St-Cyr 2007 (NCAR)

n = max(size(xgl));

% Elemental Mass matrix 1D

Mhat1D = [2/3 1/3;1/3 2/3];

% Elemental Stiffness matrix 1D

Khat1D = [1/2 -1/2;-1/2 1/2];

M = zeros(n,n);
K = zeros(n,n);

% Assembly

for i=1:n-1
    h = xgl(i+1) -xgl(i);
    M(i:i+1,i:i+1) = M(i:i+1,i:i+1) + 0.5*h*Mhat1D(:,:);
    K(i:i+1,i:i+1) = K(i:i+1,i:i+1) + (2/h)*Khat1D(:,:);
end 

Operators.M1DQ1 = M;
Operators.K1DQ1 = K;
