% Stitch schemes together given connection matrix and BC vector.
%     schmHand  - function_handle to a Scheme constructor
%     order     - order of accuracy
%     schmParam - cell array of extra parameters sent to each Scheme stored as cell arrays
%     blocks    - block definitions, On whatever form the scheme expects.
%     ms        - grid points in each direction for each block. Ex {[10,10], [10, 20]}
%     conn      - connection matrix
%     bound     - boundary condition vector, array of structs with fields w,e,s,n
%                 each field with a parameter array that is sent to schm.boundary_condition
%
% Output parameters are cell arrays and cell matrices.
%
% Ex: [schms, D, H] = stitchSchemes(schmHand, order, schmParam, blocks, ms, conn, bound)
function [schms, D, H,penalty] = stitchSchemes(schmHand, order, schmParam, grids, conn, bound,timeDep)
    default_arg('schmParam',[]);
    default_arg('timeDep','N');
    n_blocks = numel(grids);

    % Creating Schemes
    for i = 1:n_blocks
        if isempty(schmParam);
            schms{i} = schmHand(grids{i},order,[]);
        elseif ~iscell(schmParam)
            param = schmParam(i);
            schms{i} = schmHand(grids{i},order,param);
        else
            param = schmParam{i};
            if iscell(param)
                schms{i} = schmHand(grids{i},order,param{:});
            else
                schms{i} = schmHand(grids{i},order,param);
            end
        end

        % class(schmParam)
        % class(ms)
        % class(blocks)
        % class(schmParam{i})
        % class(ms)


    end


    % Total norm
    H = cell(n_blocks,n_blocks);
    for i = 1:n_blocks
        H{i,i} = schms{i}.H;
    end
    
    %% Total system matrix
    
    % Differentiation terms
    D = cell(n_blocks,n_blocks);
    penalty = cell(n_blocks,n_blocks);
    
    
    for i = 1:n_blocks
        penalty{i,i}= @(t)0;
        D{i,i} = schms{i}.D;
    end
    
    % Boundary penalty terms
    switch timeDep
        case {'n','no','N','No'}
            for i = 1:n_blocks
                if ~isstruct(bound{i})
                    continue
                end
                
                fn = fieldnames(bound{i});
                for j = 1:length(fn)
                    bc = bound{i}.(fn{j});
                    if isempty(bc)
                        continue
                    end
                    
                    [closure, ~] = schms{i}.boundary_condition(fn{j},bc{:});
                    D{i,i} = D{i,i}+closure;
                end
            end
            
            % Interface penalty terms
            for i = 1:n_blocks
                for j = 1:n_blocks
                    intf = conn{i,j};
                    if isempty(intf)
                        continue
                    end
                    
                    [uu,uv,vv,vu] = schms{i}.interface_coupling(schms{i},intf{1},schms{j},intf{2});
                    D{i,i} = D{i,i} + uu;
                    D{i,j} = uv;
                    D{j,j} = D{j,j} + vv;
                    D{j,i} = vu;
                end
            end
        case {'y','yes','Y','Yes'}
            for i = 1:n_blocks
                if ~isstruct(bound{i})
                    continue
                end
                
                fn = fieldnames(bound{i});
                for j = 1:length(fn)
                    bc = bound{i}.(fn{j});
                    if isempty(bc)
                        continue
                    end
                    
                    [closure, penaltyi] = schms{i}.boundary_condition(fn{j},bc{:});
                    D{i,i} =@(t) D{i,i}(t) + closure(t);
                    penalty{i,i} = @(t)  penalty{i,i}(t) + penaltyi(t);
                end
            end
            
            % Interface penalty terms
            for i = 1:n_blocks
                for j = 1:n_blocks
                    intf = conn{i,j};
                    if isempty(intf)
                        continue
                    end
                    
                    [uu,uv,vv,vu] = schms{i}.interface_coupling(schms{i},intf{1},schms{j},intf{2});
                    D{i,i} = @(t) D{i,i}(t) + uu(t);
                    D{i,j} = uv;
                    D{j,j} = @(t)D{j,j}(t) + vv(t);
                    D{j,i} = vu;
                end
            end
    end