function femmodel=solver_nonlinear(femmodel,conserve_loads)
%SOLVER_NONLINEAR - core solver of diagnostic run
%
%   Usage:
%      [femmodel]=solver_nonlinear(femmodel,conserve_loads)

	%Branch on partitioning schema requested
	rift_penalty_threshold=femmodel.parameters.DiagnosticRiftPenaltyThreshold;
	maxiter=femmodel.parameters.DiagnosticMaxiter;
	configuration_type=femmodel.parameters.ConfigurationType;
	[femmodel.nodes]=UpdateConstraints(femmodel.nodes,femmodel.constraints,femmodel.parameters);

	%keep a copy of loads for now
	loads=femmodel.loads;

	%initialize solution vector
	converged=0; count=1;

	%Start non-linear iteration using input velocity: 
	ug=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, loads, femmodel.materials, femmodel.parameters);
	uf=Reducevectorgtof( ug, femmodel.nodes,femmodel.parameters);

	%Update the solution to make sure that vx and vxold are similar
	[femmodel.elements loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,double(converged),ConvergedEnum);
	[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug);

	while(~converged),

		%save pointer to old velocity
		old_ug=ug;
		old_uf=uf;

		[K_ff,K_fs,p_f,df,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
		ys=CreateNodalConstraints(femmodel.nodes,configuration_type);
		p_f = Reduceload( p_f, K_fs, ys);

		issmprintf(VerboseSolver(),'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
		uf=Solver(K_ff,p_f,old_uf,df,femmodel.parameters);
		ug= Mergesolutionfromftog( uf, ys, femmodel.nodes,femmodel.parameters); 

		%Figure out if convergence have been reached
		converged=convergence(K_ff,p_f,uf,old_uf,femmodel.parameters);

		%add convergence status into  status
		[femmodel.elements loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,double(converged),ConvergedEnum);

		[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug);
		[loads,constraints_converged,num_unstable_constraints] =ConstraintsState( femmodel.elements,femmodel.nodes,femmodel.vertices,loads, femmodel.materials,femmodel.parameters);

		issmprintf(VerboseConvergence(),'%s%i','      number of unstable constraints: ',num_unstable_constraints);

		%rift convergence
		if ~constraints_converged,
			if converged,
				if num_unstable_constraints <= rift_penalty_threshold,
					converged=1;
				else 
					converged=0;
				end
			end
		end

		%increase count
		count=count+1;
		if(converged==1) break; end
		if(count>maxiter),
			issmprintf(true,'%s%i%s','      maximum number of iterations ',maxiter,' exceeded');
			break;
		end

	end

	%deal with loads:
	if conserve_loads==false,
		femmodel.loads=loads;
	end
end
