function varargout = femlib_stressWedge(varargin)
% function varargout = femlib_stressWedge(varargin)
%
%
% Dr. A. I. Hanna (2006)
error(nargchk(0,inf,nargin));
if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    error(['Arguments to ', mfilename, ' must come param/value in pairs.']);
end
Xe = [];
strain = [];
scalingMat = 1;
dof = 3;
N = 6;
GQP = [];
residualStress = [];
for i=1:2:length(varargin)
    switch lower(varargin{i})
        case 'pts'
            Xe = varargin{i+1};
        case 'gqp'
            GQP = varargin{i+1};
        case 'dof'
            dof = varargin{i+1};
        case 'n'
            N = varargin{i+1};
        case 'wgauss'
            wGauss = varargin{i+1};
        case 'materialstiffmat'
            D = varargin{i+1};
        case 'strain'
            strain = varargin{i+1};
        case 'scalingmatrix'
            scalingMat = varargin{i+1};
        otherwise
            error(['Unknown parameter name passed to ', mfilename, '.  Name was ' varargin{i}])
    end
end
if isempty(Xe)
    error('pts parameter is required.')
end
if isempty(GQP)
    error('Gaussian quadrature points parameter is required.')
end
if isempty(wGauss)
    wGauss = ones(N,1)/N;
end

if isempty(D)
    D = femlib_getMaterialStiffnessMatrix(3000, 0.35);
end
if isempty(strain)
    error('strain parameter is required.')
end
Ke = zeros(dof*N,dof*N);
fe = zeros(dof*N,1);
for i=1:N
    isoG = GQP(i,:);
    Nprime = femlib_NGrad(isoG(1), isoG(2), isoG(3));
    J = Nprime*Xe;
    G = inv(J)*Nprime;
    B = femlib_constructB(G);
    BtD = B'*D;
    if scalingMat
        scale = det(J)*wGauss(i);
    else
        scale = wGauss(i);
    end
    Ke = Ke + BtD*B*scale;
    fe = fe + BtD*strain*scale;
end
varargout{1} = Ke;
varargout{2} = fe;
return;
