%% 
%Last modified: March 18, 2014
% Start by opening the parallel toolbox

% matlabpool open     % COMMENT IF NO PARALLEL TOOLBOX

% Close at the end

% Script Experiments with HDG3d for Reaction-Diffusion problems

clear 

k=input('Polynomial degree (<=3): ');
example=input('Input example (1->P1, 2->P2, 3->P3, 4->Cinfty): ');
ctmass=input('Coefficients: (1) variable; (0) constant: ');
domain=input('Domain: (1) chimney (2) cube (3) corner: ');
switch domain
    case 1
        load('../meshes/4Tchimney')     % Chimney-shaped domain - Dirichlet = horizontal
    case 2
        load('../meshes/sixT3dDir')     % (0,1)^3 all Dirichlet
    case 3
        load('../meshes/Corner')        % Fichera corner (uniform partition) - mixex BC
end
listT={T1,T2,T3,T4};

switch ctmass
    case 1
        kappa = @(x,y,z) 2+sin(x).*sin(y).*sin(z);
        kx = @(x,y,z) cos(x).*sin(y).*sin(z);
        ky = @(x,y,z) sin(x).*cos(y).*sin(z);
        kz = @(x,y,z) sin(x).*sin(y).*cos(z);
        c = @(x,y,z) 1.+0.5*(x.^2+y.^2+z.^2);
    case 0
        kappa = @(x,y,z) 2+0.*x;
        kx = @(x,y,z) 0.*x;
        ky = @(x,y,z) 0.*x;
        kz = @(x,y,z) 0.*x;
        c = @(x,y,z) 1+0.*x;
end

formulas=checkQuadrature3d(k,ctmass);

switch example
    case 1
        u = @(x,y,z) 2*x+3*y+z;
        ux = @(x,y,z) 2+0.*x;
        uy = @(x,y,z) 3+0.*x;
        uz = @(x,y,z) 1+0.*x;
        uxx = @(x,y,z) 0.*x;
        uyy = @(x,y,z) 0.*x;
        uzz = @(x,y,z) 0.*x;
    case 2
        u = @(x,y,z) x.^2+y.*x+z.*y+2*z.^2;
        ux = @(x,y,z) 2*x + y;
        uy = @(x,y,z) x + z;
        uz = @(x,y,z) y + 4*z;
        uxx = @(x,y,z) 2+0.*x;
        uyy = @(x,y,z) 0.*x;
        uzz = @(x,y,z) 4+0.*x;
    case 3
        u = @(x,y,z) 2*x.^2.*y+4*y.^2.*z+3*x.*z.^2;
        ux = @(x,y,z) 3*z.^2 + 4*x.*y;
        uy = @(x,y,z) 2*x.^2 + 8*y.*z;
        uz = @(x,y,z) 4*y.^2 + 6*x.*z;
        uxx = @(x,y,z) 4*y;
        uyy = @(x,y,z) 8*z;
        uzz = @(x,y,z) 6*x;
    case 4
        u = @(x,y,z) sin(x.*y.*z);
        ux = @(x,y,z) y.*z.*cos(x.*y.*z);
        uy = @(x,y,z) x.*z.*cos(x.*y.*z);
        uz = @(x,y,z) x.*y.*cos(x.*y.*z);
        uxx = @(x,y,z) -y.^2.*z.^2.*sin(x.*y.*z);
        uyy = @(x,y,z) -x.^2.*z.^2.*sin(x.*y.*z);
        uzz = @(x,y,z) -x.^2.*y.^2.*sin(x.*y.*z);
end

uD = u;
km = @(x,y,z) 1./kappa(x,y,z);
qx = @(x,y,z) -kappa(x,y,z).*ux(x,y,z);
qy = @(x,y,z) -kappa(x,y,z).*uy(x,y,z);
qz = @(x,y,z) -kappa(x,y,z).*uz(x,y,z);

f = @(x,y,z) -(kx(x,y,z).*ux(x,y,z)+kappa(x,y,z).*uxx(x,y,z))...
             -(ky(x,y,z).*uy(x,y,z)+kappa(x,y,z).*uyy(x,y,z))...
             -(kz(x,y,z).*uz(x,y,z)+kappa(x,y,z).*uzz(x,y,z))...
             +c(x,y,z).*u(x,y,z);

gx = @(x,y,z) -qx(x,y,z);         
gy = @(x,y,z) -qy(x,y,z);  
gz = @(x,y,z) -qz(x,y,z); 

a = @(x,y,z) 0.*x;

ErrorU=[];
ErrorQ=[];
ErrorUhat=[];
ErrorPu=[];
ErrorPuhat=[];
ErrorUstar=[];
h=[];

% Norms of unknowns for relative error

Tmax=listT{end};
Nelts=size(Tmax.elements,1);
d3=nchoosek(k+3,3);
normU=errorElem(Tmax,u,zeros(d3,Nelts),k,formulas{1});
normQ=errorElem(Tmax,qx,zeros(d3,Nelts),k,formulas{1})...
        +errorElem(Tmax,qy,zeros(d3,Nelts),k,formulas{1})...
        +errorElem(Tmax,qz,zeros(d3,Nelts),k,formulas{1});
    
for i=1:length(listT)
    T=listT{i};
    Nelts=size(T.elements,1);
    tau=createTau3d(Nelts,2);    
    [Uh,Qxh,Qyh,Qzh,Uhat]=HDG3d(km,c,f,tau,T,k,formulas,uD,gx,gy,gz);
    normUhat=errorFaces(T,u,zeros(size(Uhat)),k,formulas{4});
    
    % Errors with exact solution
    error_uhat=errorFaces(T,u,Uhat,k,formulas{4});
    error_q   =errorElem(T,qx,Qxh,k,formulas{1})...
               +errorElem(T,qy,Qyh,k,formulas{1})...
               +errorElem(T,qz,Qzh,k,formulas{1});
    error_u   =errorElem(T,u,Uh,k,formulas{1});
    ErrorQ=[ErrorQ error_q/normQ]; 
    ErrorU=[ErrorU error_u/normU];
    ErrorUhat=[ErrorUhat error_uhat/normUhat];
    
    % Errors with projections
    Puhat=L2projskeleton3d(u,T,k,formulas{3},formulas{4});
    [Pqx,Pqy,Pqz,Pu]=projectHDG3d(T,{qx,qy,qz,u},k,tau,formulas);
    error_Puhat=errorFaces(T,a,Uhat-Puhat,k,formulas{4});
    error_Pu=errorElem(T,a,Pu-Uh,k,formulas{1});
    ErrorPu=[ErrorPu error_Pu/normU];
    ErrorPuhat=[ErrorPuhat error_Puhat/normUhat];
    
    % Postprocessing
    Uhstar=postprocessing(T,km,Qxh,Qyh,Qzh,Uh,k,formulas{1});
    error_Ustar=errorElem(T,u,Uhstar,k+1,formulas{1});
    ErrorUstar=[ErrorUstar error_Ustar/normU];
       
    h=[h 1/(2^i)];
end

rateQ=log2(ErrorQ(1:end-1)./ErrorQ(2:end));
rateU=log2(ErrorU(1:end-1)./ErrorU(2:end));
rateUhat=log2(ErrorUhat(1:end-1)./ErrorUhat(2:end));
ratePu=log2(ErrorPu(1:end-1)./ErrorPu(2:end));
ratePuhat=log2(ErrorPuhat(1:end-1)./ErrorPuhat(2:end));
rateUstar=log2(ErrorUstar(1:end-1)./ErrorUstar(2:end));

format bank
disp('Rate |Q-Qh|   Rate |U-Uh|     Rate |U-Uhat|_h');
disp([rateQ' rateU' rateUhat']);
disp('Rate |\Pi u-Uh|  Rate |Pu-Uhat|_h  Rate |U-U*|')
disp([ratePu' ratePuhat' rateUstar']);

format shortE
disp('Error |Q-Qh|   Error |U-Uh|     Error |U-Uhat|_h');
disp([ErrorQ' ErrorU' ErrorUhat']);
disp('Error |\Pi u-Uh|  Error |Pu-Uhat|_h Error |U-U*|')
disp([ErrorPu' ErrorPuhat' ErrorUstar']);





%%

% matlabpool close        % COMMENT IF NO PARALLEL TOOLBOX