function melting_g=icemelting_core(m,analysis_type,varargin);
%ICEMELTING_CORE - core solution of melting computation
%
%   Usage:
%      melting_g=icemelting_core(m,analysis_type,varargin)

%global variables
global cluster gridset

%recover fem model fields
elements=m.elements;
grids=m.grids;
materials=m.materials;
loads=m.loads;
ys=m.ys;
gridset=m.gridset;
G_mn=m.Gmn;
params=m.params;

%recover parameters
sparsity=params.sparsity;
solver_type=params.solver_type;

%recover extra inputs and initialize solution
if nargin==3,
	inputs=varargin{1};
else	
	inputs=struct();
end

%initialize solution vector
melting_g=[];

disp(sprintf('%s','   solving'));

% Generate system matrices (stiffness and load)
%stiffness is computed at each step, not the loads.
kflag=1;pflag=1;

%generate stiffness and loads
[K_gg,p_g]=SystemMatrices(elements,grids,loads,materials,kflag, pflag, params.sparsity,inputs,'melting');
[K_gg,p_g]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, params.sparsity,inputs,'melting',K_gg,p_g);

% Reduce stiffness matrix from g set to f set (effectively rooting out the single point constraints),
% and compute modifications to loads from single point constraints.
[K_ff K_fs]=IceReducematrixfromgtof(K_gg,G_mn); 

% Reduce load from g set to f set
p_f=IceReducerightside(p_g,G_mn,K_fs,ys);

% Solve
melting_f=Solver(K_ff,p_f,params.solver_type);

% Add single point constraints back, ie increase f-set by s-set into the global g-set.
melting_g= Mergesolution_g(melting_f,G_mn,ys); 
