Index: /issm/trunk/src/m/solutions/NewFemModel.m
===================================================================
--- /issm/trunk/src/m/solutions/NewFemModel.m	(revision 5775)
+++ /issm/trunk/src/m/solutions/NewFemModel.m	(revision 5776)
@@ -19,5 +19,4 @@
    
 	%Initialize some fiels with empty celils
-	femmodel.m_yg=cell(nummodels,1);
 	femmodel.m_nodesets=cell(nummodels,1);
 	femmodel.m_ys=cell(nummodels,1);
@@ -32,20 +31,22 @@
 		femmodel=SetCurrentConfiguration(femmodel,analysis_type);
 
-		displaystring(md.verbose,'%s','      generating degrees of freedom...');
+		displaystring(md.verbose,'%s','      generating vertices degrees of freedom');
 		if ~isfield(femmodel,'part'),
 			[femmodel.vertices,femmodel.part,femmodel.tpart]=VerticesDof(femmodel.vertices, femmodel.parameters); %do not create partition vector twice! we only have one set of vertices!
 		end
+
+		displaystring(md.verbose,'%s','      resolve node constraints');
+		[femmodel.nodes]=SpcNodes(femmodel.nodes,femmodel.constraints,analysis_type);
+
+		displaystring(md.verbose,'%s','      create nodal degrees of freedom');
 		[femmodel.nodes]=NodesDof(femmodel.nodes,femmodel.parameters);
 
-		displaystring(md.verbose,'%s','      generating single point constraints...');
-		[femmodel.nodes,femmodel.m_yg{i}]=SpcNodes(femmodel.nodes,femmodel.constraints,analysis_type);
+		displaystring(md.verbose,'%s','      create nodal constraints vector');
+		femmodel.m_ys{i}=CreateNodalConstraints(femmodel.nodes,analysis_type);
 
-		displaystring(md.verbose,'%s','      generating node sets...');
+		displaystring(md.verbose,'%s','      create node sets');
 		femmodel.m_nodesets{i}=BuildNodeSets(femmodel.nodes,analysis_type);
 
-		displaystring(md.verbose,'%s','      reducing single point constraints vector...');
-		femmodel.m_ys{i}=Reducevectorgtos(femmodel.m_yg{i},femmodel.m_nodesets{i});
-
-		displaystring(md.verbose,'%s','      configuring elements and loads...');
+		displaystring(md.verbose,'%s','      configuring elements and loads');
 		[femmodel.elements,femmodel.loads,femmodel.nodes,femmodel.parameters] = ConfigureObjects( femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices,femmodel.materials,femmodel.parameters);
 	end
Index: /issm/trunk/src/m/solutions/ResetBoundaryConditions.m
===================================================================
--- /issm/trunk/src/m/solutions/ResetBoundaryConditions.m	(revision 5775)
+++ /issm/trunk/src/m/solutions/ResetBoundaryConditions.m	(revision 5776)
@@ -20,13 +20,10 @@
 	%For this analysis_type, free existing boundary condition vectors:
 	analysis_counter=femmodel.parameters.AnalysisCounter+1; %matlab indexing on counter
-
-	%global dof set
-	femmodel.m_yg{analysis_counter}=[];
-	%in the s-set
 	femmodel.m_ys{analysis_counter}=[];
 
-	%Now, duplicate ug (the solution vector) into the boundary conditions vector on the g-set
-	femmodel.m_yg{analysis_counter}=ug;
+	%Reduce from g to s set
+	ys=Reducevectorgtos(ug,femmodel.m_nodesets{analysis_counter});
 
-	%Reduce from g to s set
-	femmodel.m_ys{analysis_counter}=Reducevectorgtos(femmodel.m_yg{analysis_counter},femmodel.m_nodesets{analysis_counter});
+	%in the s-set
+	femmodel.m_ys{analysis_counter}=ys;
+
Index: /issm/trunk/src/m/solutions/SetCurrentConfiguration.m
===================================================================
--- /issm/trunk/src/m/solutions/SetCurrentConfiguration.m	(revision 5775)
+++ /issm/trunk/src/m/solutions/SetCurrentConfiguration.m	(revision 5776)
@@ -39,5 +39,4 @@
 	%activate matrices and vectors: 
 	femmodel.nodesets=femmodel.m_nodesets{found};
-	femmodel.yg=femmodel.m_yg{found};
 	femmodel.ys=femmodel.m_ys{found};
 
Index: /issm/trunk/src/m/solvers/solver_adjoint_linear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_adjoint_linear.m	(revision 5775)
+++ /issm/trunk/src/m/solvers/solver_adjoint_linear.m	(revision 5776)
@@ -5,5 +5,5 @@
 %      femmodel =solver_adjoint_linear(femmodel)
 
-	[K_gg, p_g, kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+	[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 	
 	[K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters); 
Index: /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m	(revision 5775)
+++ /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m	(revision 5776)
@@ -27,5 +27,5 @@
 		old_uf=uf;
 		
-		[K_gg, p_g, kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
+		[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
 		
 		[K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters); 
Index: /issm/trunk/src/m/solvers/solver_linear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_linear.m	(revision 5775)
+++ /issm/trunk/src/m/solvers/solver_linear.m	(revision 5776)
@@ -5,5 +5,5 @@
 %      femmodel =solver_linear(femmodel)
 
-	[K_gg, p_g,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+	[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 	
 	[K_ff, K_fs] = Reducematrixfromgtof( K_gg,  femmodel.nodesets,femmodel.parameters); 
Index: /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m	(revision 5775)
+++ /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m	(revision 5776)
@@ -20,5 +20,5 @@
 	while(~converged),
 
-		[K_gg,p_g,melting_offset]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+		[K_gg,K_ff,K_fs,p_g,p_f,melting_offset]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 
 		[K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters);