% FEm.PSE template script , multi-material, with sources, quadratic basis
global X1;        % node coordinates

% set up problem parameters
Tb = 0;           % fixed temperature applied to both ends
Aleft = 100;      % source strength in left half-domain
Aright = 10;      % source strength in right half-domain
L = 1;            % half-length of the domain
kleft = 0.1;      % conductivity in left half-domain
kright = 0.001;   % conductivity in left half-domain
nnodes = 9;       % 9 grid points for four quadratic elements
nelements = 4;    % 4 elements

% set up grid  data
X1 = linspace(0,2*L,nnodes); % grid data

% set up nodal source data
SRC = (abs(sin(pi*X1/L)))'; %  uses symmetry of sine for right half-domain
nnhd = (nnodes-1)/2; % number of nodes in a half-domain
SRC(1:nnhd,1) = SRC(1:nnhd,1)*Aleft; % apply Aleft to left half-domain
SRC(nnhd+2:nnodes,1) = SRC(nnhd+2:nnodes,1)*Aright; % apply Aright

% set up element conductivities
nel2 = nelements/2; % number of elements in each half-domain
elcond = ones(nelements,1)*kright;%create elcond with all values set to kright
elcond(1:nel2,1) = ones(nel2,1)*kleft; % put in kleft


load femlib;  % load library matrix data

% set up diffusion matrix
DIFF = asjac1D([],elcond,[],-1,'A211Q',[]);

% set up data column matrix
M = asjac1D([],[],[],1,'A200Q',[]);
b = M * SRC;

% enforce the left Dirichlet B.C.
DIFFdiri = DIFF;  % copy DIFF
DIFFdiri(1,:) = zeros(1,nnodes);
DIFFdiri(1,1) = 1;
b(1,1) = Tb;

%  enforce the right Dirichlet B.C.
DIFFdiri(nnodes,:) = zeros(1,nnodes);
DIFFdiri(nnodes,nnodes) = 1;
b(nnodes,1) = Tb;

%  Solve for Q
Q = DIFFdiri \ b

% Calculate the energy norm
format long; % display Enorm in long format
Enorm = 0.5 * Q' * DIFF * Q
format short;


