function P = ntw_cgikl(varargin)
% NTW_CGI Construct an example from Nobile et al. with decaying coefficients
%
% This function requies the PDE toolbox.  It also produces a warning when
% we generate the mesh, you can safely ignore this warning.
%
% P = ntw_cgikl(); constructs a problem with 6 parameters
% P = ntw_cgikl(d); constructs a problem with d parameters
% P = ntw_cgikl('dim',d); constructs a problem with d parameters
%
% Return:
% The output P is a struct with the following fields:
% P.A - a parameterized matrix interface
% P.Av - a parameterized matrix-vector interface
% P.b - a parameterized right hand side
% P.N - the size of the matrix
% P.d - the number of parameters (1) for this problem
% P.solve - a parameterized solution
% P.s - the vector of parameters for this problem
% P.mesh - details on the mesh (points, edges, triangles);
%
% Example:
%   P = ntw_cgikl();
%   [X,r] = spectral_galerkin(P.Av,P.b,P.s,2);
%   scatter(P.mesh.p(:,1),P.mesh.p(:,2),1,X.coefficient(:,1));
%

% Copyright, Stanford University, 2010
% Paul G. Constantine, David F. Gleich

d = [];
if nargin==1,
    d = varargin{1};
else
    args = struct(varargin{:});
    if isfield(args,'dim'), d=args.dim; dim.size = []; end
end
if isempty(d), d=6; end

% Geometry and boundary conditions
g=[2 2 2 2;0 1 1 0;1 1 0 0;0 0 1 1;0 1 1 0;1 1 1 1;0 0 0 0];
[p,e,t]=initmesh(g,'Hmax',0.03);

b=[1 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1;48 48 48 48;
    48 48 48 48;49 49 49 49;48 48 48 48];

Data = load('sisc_kl.mat'); % Loads KL
KL = Data.KL;

% set these as nested functions so they share all the loaded data
    function [K,F]=ntw_assemble(s) 
        if length(s)~=d, error('pmpack:wrongParameterSize',...
                'this ntw problem takes %i parameters, not %i', d, length(s)); end

        c=@(p,t,u,time) sisc_coeff(p,t,KL,s);

        % Solve the system.
        [K,F]=assempde(b,p,e,t,c,0,1);
    end

    function K=ntw_mat(s)
        K = ntw_assemble(s);
    end

% there is no parameter dependence on the right hand side, so just
% precompute it and make it constant.

[K0,F0] = ntw_assemble(zeros(1,d)); %#ok<ASGLU>

    function F=ntw_rhs(s) %#ok<INUSD>
        F = F0;
    end

    function Av=ntw_matvec(s,v)
        K = ntw_assemble(s);
        Av = K*v;
    end

    function r=ntw_solve(s)
        [K,F] = ntw_assemble(s);
        r = K\F;
    end

P.b = @ntw_rhs;
P.A = @ntw_mat;
P.Av = @ntw_matvec;
P.solve =@ntw_solve;
P.b_const = 1;
P.d = d;
P.N = length(p);    
variables(d) = jacobi_parameter;
for i=1:d-1
    variables(i) = jacobi_parameter;
end
P.s = variables;
P.mesh.p = p;
P.mesh.e = e;
P.mesh.t = t;

% this function is 
end

% NTW_COEFF returns a vector with the value of the NTW elliptic coefficient
% given at the cell centers of the mesh.
function ct=sisc_coeff(p,t,KL,s)

d=length(s);
s=2*KL.sv(1:d).*s(:);

cn=exp(KL.mean+KL.bases(:,1:d)*s);

ct=pdeintrp(p,t,cn);

end
