function [t_g loads melting_offset]=icethermal_core(m,analysis_type,varargin);
%ICETHERMAL_CORE - core of the thermal solution
%
%   Usage:
%      [t_g loads melting_offset]=icethermal_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
t_g=[];

converged=0;
firstiteration=1;

disp(sprintf('%s','   starting direct shooting method'));

dt=inputs.dt;

%Recover pressure from inputs
pressure=inputs.pressure;

while(~converged),
 
	% Generate system matrices (stiffness and load)
    
	%loads and stiffness are computed only once, not the penalties
	kflag=1;pflag=1;

	%generate stiffness and loads
	if firstiteration
		[K_gg,p_g]=SystemMatrices(elements,grids,loads,materials,kflag, pflag, params.sparsity,inputs,analysis_type);
		firstiteration=0;
	end
	[K_gg_penalized,p_g_penalized melting_offset]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, params.sparsity,inputs,analysis_type,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_penalized,G_mn); 

	% Reduce load from g set to f set
	p_f=IceReducerightside(p_g_penalized,G_mn,K_fs,ys);
	
	% Solve
	t_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.
	t_g= Mergesolution_g(t_f,G_mn,ys); 

	[loads,converged]=IcePenaltyConstraints(loads,grids,materials,params,struct('pressure',pressure, 'dt',dt, 'temperature',t_g));

end
