function meshInfo = femlib_wedgeIteration(meshInfo)
% function X = femlib_wedgeIteration(X, meshInfo)
%
% A script to perform a single iteration of the FEM procedure for growing a
% body defined by wedge elements.
%
% Dr. A. I. Hanna (2006)
N = length(meshInfo.wGauss);
K = zeros(meshInfo.nodes*meshInfo.dof);
f = zeros(meshInfo.nodes*meshInfo.dof,1);
if ~isfield(meshInfo, 'scalingMatrix')
    meshInfo.scalingMatrix = 1;
end
X = meshInfo.X;
Xold = X;
for i=1:length(meshInfo.element)
    ind = meshInfo.element(i).ind;
    Xe = X([ind, ind+meshInfo.dorsalDivide],:);
    strain = meshInfo.element(i).preStrain.*mean(meshInfo.vertex_temp(ind));
    [Ke, fe] = femlib_stressWedge('pts', Xe, 'GQP', meshInfo.GQP, 'dof', meshInfo.dof, 'N', N, 'wgauss', meshInfo.wGauss, 'materialstiffmat', meshInfo.D, 'Strain', strain, 'scalingMatrix', meshInfo.scalingMatrix);
    ind = [ind*3-2; ind*3-1; ind*3]; 
    ind = ind(:)';
    ind  = cat(2, ind, ind + (size(X,1)/2)*3);
    K(ind, ind) = K(ind, ind) + Ke;
    f(ind) = f(ind) - fe;
end
[K, f, usedEquations] = femlib_eliminateEquations( K, f, meshInfo.fixedEquations);
tol = 1e-4;
maxiter = size(K,1)*50;
%[u, cgflag, cgrelres, cgiter] = cgs(K, f, tol, maxiter);
[u, cgflag, cgrelres, cgiter] = qmr(K, f, tol, maxiter);


fprintf(' - CGS converged in %d iterations\n', cgiter);
if cgiter==maxiter
    fprintf(' - reached maximum iteration number in CGS\n');
end
U = zeros(meshInfo.nodes*meshInfo.dof,1);
U(usedEquations) = u;
u = reshape(U, 3, length(U)/3)';
X = X + u;
%max(abs(X(:,3)))
meshInfo.X = X ;
return;
