%%
%% Script: ML_FEM_1D_DiffusionVariableCoeff_MMS.m
%%
%% FEM code for 1D diffusion equation with variable coefficient diffusivity. 
%% Demonstration of the Method of Manufactured Solutions (MMS)
%% 26.07.2011, Dave May, ETH Zurich
%%


function [] = main()

clear all, close all, clc

% Numerical parameters
% Mesh sequence, no_elements = { 4, 8, 16, 32, 64, 128, 256, 512, 1024 }
no_el       = 128;
no_nodes    = no_el + 1;
no_nodes_el = 2;
width       = 2;
X           = [-1:width/(no_el):-1+width];
dx          = X(2)-X(1);
% Initialization of matrices and vectors
KK  = zeros(no_nodes,no_nodes);
u   = zeros(no_nodes,1);
F   = zeros(no_nodes,1);

% Node numbering
for i=1:no_el
    iel = i; 
    nodes(iel,1) = i; 
    nodes(iel,2) = i+1;
end

% Setup system matrix
for iel=1:no_el
    ileft  = nodes(iel,1);
    iright = nodes(iel,2);
    x_centroid = 0.5 * ( X(iright) + X(ileft) );

    % Physical parameters evaluated at the element centroid
    A   = MMS_Coefficient( x_centroid );
    B   = MMS_RightHandSide( x_centroid );
    
    % Setup element matrix
    K_loc   = [-A/dx  A/dx; A/dx  -A/dx];   % Matrix for conductive terms, see Maple script
    F_loc   = [B*dx/2 B*dx/2];              % Right hand side vector
    
    for i=1:no_nodes_el
        ii                              = nodes(iel,i);
        for j=1:no_nodes_el
            jj                          = nodes(iel,j);
            KK(ii,jj)                   = KK(ii,jj)       + K_loc(i,j);
        end
        F(ii)                           = F(ii)           + F_loc(i);
    end
end

% Boundary conditions
% Evaluate the manufactured solution at the end points
KK(1,:)   = 0; KK(1,1)     = 1; F(1)   = MMS_Solution( X(1)   ); % Left  boundary 
KK(end,:) = 0; KK(end,end) = 1; F(end) = MMS_Solution( X(end) ); % Right boundary

% Solve matrix
u                                       = KK\F;

% Compute the L2 error using a 1-point quadrature over each element
L2_error = 0.0;
for iel=1:no_el
    ileft  = nodes(iel,1);
    iright = nodes(iel,2);

    dV             = X(iright) - X(ileft);            % element volume
    x_centroid     = 0.5 * ( X(iright) + X(ileft) );  % coordinate of element centroid
    T_fem_centroid = 0.5 * ( u(iright) + u(ileft) );  % FE solution at centroid
    T_mms          = MMS_Solution( x_centroid );      % MMS at centroid
    T_error        = T_mms - T_fem_centroid;          % Discretisation error at centroid
    
    L2_error = L2_error + T_error * T_error * dV;    
end
L2_error = sqrt( L2_error );
fprintf(1,'h =          E_2 =  \n' );
fprintf(1,'%1.4e,  %1.6e ; \n', dx, L2_error );



% Evaluate the analytical solution and plot
L                                       = width;
x                                       = [-1:dx/2.0:-1+L];
for i=1:length(x)
    u_mms(i)  = MMS_Solution( x(i) );
end
    
% Graphics
plot(X,u,'-ok',x,u_mms,'-b'); 
xlabel('x'); ylabel('u'); grid on; 
legend('finite element soln.','manufactured soln.',4); drawnow;


% =========================================================================
%
% Manufactured coefficients, solution and right hand side
%
function val = MMS_Coefficient(x)
   val = ( 4.0 * tanh(10.0 * x) + 6.0 ) * exp(0.8 * x);

function val = MMS_Solution(x)
   val = 8.0 * sin(x)*x^4  + 20.0;

function val = MMS_RightHandSide(x)
   val = (6 + 4*tanh(10*x))*(-8*x^4*sin(x) + 64*x^3*cos(x)                ...
   + 96*x^2*sin(x))*exp(0.8*x) + (40 - 40*tanh(10*x)^2)*(8*x^4*cos(x)     ...
   + 32*x^3*sin(x))*exp(0.8*x) + 0.8*(6 + 4*tanh(10*x))*(8*x^4*cos(x)     ...
   + 32*x^3*sin(x))*exp(0.8*x);

   
