Index: /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26681)
+++ /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26682)
@@ -136,8 +136,5 @@
 	end
 
-	print(pe)
-	error("STOP")
-
-	return Pe
+	return pe
 end #}}}
 function GetSolutionFromInputs(analysis::StressbalanceAnalysis,ug::Vector{Float64},element::Tria) #{{{
Index: /issm/trunk-jpl/src/jl/solve/elementmatrix.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/elementmatrix.jl	(revision 26681)
+++ /issm/trunk-jpl/src/jl/solve/elementmatrix.jl	(revision 26682)
@@ -31,4 +31,35 @@
 	display(this.values)
 end# }}}
+function AddToGlobal!(Ke::ElementMatrix,Kff::IssmMatrix,Kfs::IssmMatrix)#{{{
+
+	#First check that the element matrix looks alright
+	CheckConsistency(Ke)
+
+	#See if we need to do anything
+	is_fset = false
+	is_sset = false
+	for i in 1:Ke.nrows
+		if(Ke.fglobaldoflist[i]>0) is_fset = true end
+		if(Ke.sglobaldoflist[i]>0) is_sset = true end
+	end
+
+	if is_fset
+		AddValues!(Kff,Ke.nrows,Ke.fglobaldoflist,Ke.nrows,Ke.fglobaldoflist,Ke.values)
+	end
+	if is_sset
+		AddValues!(Kfs,Ke.nrows,Ke.fglobaldoflist,Ke.nrows,Ke.sglobaldoflist,Ke.values)
+	end
+
+end#}}}
+function CheckConsistency(Ke::ElementMatrix)#{{{
+
+	for i in 1:Ke.nrows
+		for j in 1:Ke.nrows
+			if(isnan(Ke.values[i,j])) error("NaN found in Element Matrix") end
+			if(isinf(Ke.values[i,j])) error("Inf found in Element Matrix") end
+			if(abs(Ke.values[i,j])>1.e+50) error("Element Matrix values exceeds 1.e+50") end
+		end
+	end
+end#}}}
 
 mutable struct ElementVector#{{{
@@ -58,2 +89,29 @@
 	display(this.values)
 end# }}}
+function AddToGlobal!(pe::ElementVector,pf::IssmVector)#{{{
+
+	#First check that the element matrix looks alright
+	CheckConsistency(pe)
+
+	#See if we need to do anything
+	is_fset = false
+	for i in 1:pe.nrows
+		if(pe.fglobaldoflist[i]>0)
+			is_fset = true 
+			break
+		end
+	end
+
+	if is_fset
+		AddValues!(pf,pe.nrows,pe.fglobaldoflist,pe.values)
+	end
+
+end#}}}
+function CheckConsistency(pe::ElementVector)#{{{
+
+	for i in 1:pe.nrows
+		if(isnan(pe.values[i])) error("NaN found in Element Vector") end
+		if(isinf(pe.values[i])) error("Inf found in Element Vector") end
+		if(abs(pe.values[i])>1.e+50) error("Element Vector values exceeds 1.e+50") end
+	end
+end#}}}
Index: /issm/trunk-jpl/src/jl/solve/modules.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/modules.jl	(revision 26681)
+++ /issm/trunk-jpl/src/jl/solve/modules.jl	(revision 26682)
@@ -171,13 +171,25 @@
 
 	#Construct Stiffness matrix and load vector from elements
-	for i=1:length(femmodel.elements)
+	for i in 1:length(femmodel.elements)
 		Ke = CreateKMatrix(analysis,femmodel.elements[i])
 		pe = CreatePVector(analysis,femmodel.elements[i])
-		print(Ke)
-		error("don't know what to do")
+
+		AddToGlobal!(Ke,Kff,Kfs)
+		AddToGlobal!(pe,pf)
 	end
 	
-	
-	
-	error("STOP")
+	return Kff, Kfs, pf
 end# }}}
+function CreateNodalConstraintsx(nodes::Vector{Node})# {{{
+
+	#Allocate vector
+	ssize=NumberOfDofs(nodes,SsetEnum)
+	ys=IssmVector(ssize)
+
+	#constraints vector with the constraint values
+	for i in 1:length(nodes)
+		CreateNodalConstraints(nodes[i],ys)
+	end
+
+	return ys
+end# }}}
Index: /issm/trunk-jpl/src/jl/solve/nodes.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/nodes.jl	(revision 26681)
+++ /issm/trunk-jpl/src/jl/solve/nodes.jl	(revision 26682)
@@ -30,4 +30,11 @@
 	node.sdoflist[dof]  = +1
 	node.svalues[dof]   = value
+
+end# }}}
+function CreateNodalConstraints(node::Node,ys::IssmVector) #{{{
+
+	if(SSize(node)>0)
+		SetValues!(ys,node.gsize,node.sdoflist,node.svalues)
+	end
 
 end# }}}
@@ -138,4 +145,17 @@
 
 end# }}}
+function SSize(node::Node) #{{{
+
+	ssize = 0
+
+	for i=1:node.gsize
+		if node.fdoflist[i]!=-1
+			ssize+=1
+		end
+	end
+
+	return ssize
+
+end# }}}
 
 #Nodes functions
Index: /issm/trunk-jpl/src/jl/solve/solutionsequences.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/solutionsequences.jl	(revision 26681)
+++ /issm/trunk-jpl/src/jl/solve/solutionsequences.jl	(revision 26682)
@@ -17,4 +17,8 @@
 		#Get new matrices
 		Kff, Kfs, pf = SystemMatricesx(femmodel,analysis)
+
+		#Enforce constraints
+		ys = CreateNodalConstraintsx(femmodel.nodes)
+
 		error("not implemented yet")
 
Index: /issm/trunk-jpl/src/jl/solve/solve.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/solve.jl	(revision 26681)
+++ /issm/trunk-jpl/src/jl/solve/solve.jl	(revision 26682)
@@ -1,4 +1,5 @@
 include("../md/classes.jl")
 include("./issmenums.jl")
+include("./toolkits.jl")
 include("./gauss.jl")
 include("./parameters.jl")
@@ -14,7 +15,6 @@
 include("./solutionsequences.jl")
 include("./modules.jl")
-include("./toolkits.jl")
+include("./elementmatrix.jl")
 include("./utils.jl")
-include("./elementmatrix.jl")
 
 function IssmCore(md::model) #{{{
Index: /issm/trunk-jpl/src/jl/solve/toolkits.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/toolkits.jl	(revision 26681)
+++ /issm/trunk-jpl/src/jl/solve/toolkits.jl	(revision 26682)
@@ -8,4 +8,17 @@
 	return IssmMatrix(spzeros(M,N))
 end#}}}
+function AddValues!(matrix::IssmMatrix,m::Int64,midx::Vector{Int64},n::Int64,nidx::Vector{Int64},values::Matrix{Float64})#{{{
+
+	#This is inefficient now, but it will work
+	for i in 1:m
+		if(midx[i]==-1) continue end
+		for j in 1:n
+			if(nidx[j]==-1) continue end
+			matrix.matrix[midx[i],nidx[j]] += values[i,j]
+		end
+	end
+
+end#}}}
+
 mutable struct IssmVector #{{{
 	vector::Vector{Float64}
@@ -14,2 +27,20 @@
 	return IssmVector(zeros(M))
 end#}}}
+function AddValues!(vector::IssmVector,m::Int64,midx::Vector{Int64},values::Vector{Float64})#{{{
+
+	#This is inefficient now, but it will work
+	for i in 1:m
+		if(midx[i]==-1) continue end
+		vector.vector[midx[i]] += values[i]
+	end
+
+end#}}}
+function SetValues!(vector::IssmVector,m::Int64,midx::Vector{Int64},values::Vector{Float64})#{{{
+
+	#This is inefficient now, but it will work
+	for i in 1:m
+		if(midx[i]==-1) continue end
+		vector.vector[midx[i]] = values[i]
+	end
+
+end#}}}