Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 13876)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 13877)
@@ -299,6 +299,4 @@
 					./modules/StringToEnumx/StringToEnumx.cpp\
 					./modules/StringToEnumx/StringToEnumx.h\
-					./modules/SystemMatricesx/SystemMatricesx.cpp\
-					./modules/SystemMatricesx/SystemMatricesx.h\
 					./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp\
 					./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h\
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 13876)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 13877)
@@ -518,4 +518,108 @@
 }
 /*}}}*/
+void FemModel::SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax){/*{{{*/
+
+	/*intermediary: */
+	int      i;
+	int      fsize,ssize;
+	int      connectivity, numberofdofspernode;
+	int      analysis_type, configuration_type;
+	Element *element = NULL;
+	Load    *load    = NULL;
+
+	/*output: */
+	Matrix<IssmDouble> *Kff  = NULL;
+	Matrix<IssmDouble> *Kfs  = NULL;
+	Vector<IssmDouble> *pf   = NULL;
+	Vector<IssmDouble> *df   = NULL;
+	IssmDouble          kmax;
+
+	/*Display message*/
+	if(VerboseModule()) _pprintLine_("   Generating matrices");
+
+	/*retrive parameters: */
+	this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+	this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	this->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
+
+	/*Get size of matrices: */
+	fsize=this->nodes->NumberOfDofs(configuration_type,FsetEnum);
+	ssize=this->nodes->NumberOfDofs(configuration_type,SsetEnum);
+
+	numberofdofspernode=this->nodes->MaxNumDofs(configuration_type,GsetEnum);
+
+	/*Compute penalty free mstiffness matrix and load vector*/
+	Kff=new Matrix<IssmDouble>(fsize,fsize,connectivity,numberofdofspernode);
+	Kfs=new Matrix<IssmDouble>(fsize,ssize,connectivity,numberofdofspernode);
+	df =new Vector<IssmDouble>(fsize);
+
+	/*Fill stiffness matrix from elements: */
+	for (i=0;i<this->elements->Size();i++){
+		element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		element->CreateKMatrix(Kff,Kfs,df);
+	}
+
+	/*Fill stiffness matrix from loads if loads have the current configuration_type: */
+	for (i=0;i<this->loads->Size();i++){
+		load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
+		if (load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs);
+	}
+
+	/*Assemble matrices*/
+	Kff->Assemble();
+	Kfs->Assemble();
+	df->Assemble();
+
+	/*Create Load vector*/
+	pf=new Vector<IssmDouble>(fsize);
+
+	/*Fill right hand side vector, from elements: */
+	for (i=0;i<this->elements->Size();i++){
+		element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		element->CreatePVector(pf);
+	}
+
+	/*Fill right hand side from loads if loads have the current configuration_type: */
+	for (i=0;i<this->loads->Size();i++){
+		load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
+		if (load->InAnalysis(configuration_type)) load->CreatePVector(pf);
+	}
+	pf->Assemble();
+
+	/*Now, figure out maximum value of K_gg, so that we can penalize it correctly: */
+	kmax=Kff->Norm(NORM_INF);
+
+	/*Now, deal with penalties*/
+	/*Fill stiffness matrix from loads: */
+	for (i=0;i<this->loads->Size();i++){
+		load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
+		if (load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
+	}
+
+	/*Assemble matrices*/
+	Kff->Assemble();
+	Kfs->Assemble();
+
+	/*Fill right hand side vector, from loads: */
+	for (i=0;i<this->loads->Size();i++){
+		load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
+		if (load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pf,kmax);
+	}
+
+	/*Assemble vector*/
+	pf->Assemble();
+
+	/*Assign output pointers: */
+	if(pKff) *pKff=Kff;
+	else      xdelete(&Kff);
+	if(pKfs) *pKfs=Kfs;
+	else      xdelete(&Kfs);
+	if(ppf)  *ppf=pf;
+	else      xdelete(&pf);
+	if(pdf)  *pdf=df;
+	else      xdelete(&df);
+	if(pkmax) *pkmax=kmax;
+}
+/*}}}*/
 void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 13876)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 13877)
@@ -89,4 +89,5 @@
 		void ThicknessAbsGradientx( IssmDouble* pJ, bool process_units,int weight_index);
 		#endif
+		void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax);
 		void TimeAdaptx(IssmDouble* pdt);
 		void UpdateConstraintsx(void);
Index: /issm/trunk-jpl/src/c/modules/modules.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/modules.h	(revision 13876)
+++ /issm/trunk-jpl/src/c/modules/modules.h	(revision 13877)
@@ -91,5 +91,4 @@
 #include "./SpcNodesx/SpcNodesx.h"
 #include "./SurfaceAreax/SurfaceAreax.h"
-#include "./SystemMatricesx/SystemMatricesx.h"
 #include "./CreateJacobianMatrixx/CreateJacobianMatrixx.h"
 #include "./TriaSearchx/TriaSearchx.h"
Index: /issm/trunk-jpl/src/c/solvers/solver_adjoint_linear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_adjoint_linear.cpp	(revision 13876)
+++ /issm/trunk-jpl/src/c/solvers/solver_adjoint_linear.cpp	(revision 13877)
@@ -26,5 +26,5 @@
 	femmodel->UpdateConstraintsx();
 
-	SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+	femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
 	CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 	Reduceloadx(pf, Kfs, ys,true); xdelete(&Kfs); //true means spc = 0
Index: /issm/trunk-jpl/src/c/solvers/solver_linear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_linear.cpp	(revision 13876)
+++ /issm/trunk-jpl/src/c/solvers/solver_linear.cpp	(revision 13877)
@@ -24,5 +24,5 @@
 	femmodel->UpdateConstraintsx();
 
-	SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+	femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
 	CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 	Reduceloadx(pf, Kfs, ys); xdelete(&Kfs);
Index: /issm/trunk-jpl/src/c/solvers/solver_newton.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 13876)
+++ /issm/trunk-jpl/src/c/solvers/solver_newton.cpp	(revision 13877)
@@ -56,5 +56,5 @@
 
 		/*Solver forward model*/
-		SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		Reduceloadx(pf,Kfs,ys);xdelete(&Kfs);
@@ -83,5 +83,5 @@
 
 		/*Prepare next iteration using Newton's method*/
-		SystemMatricesx(&Kff,&Kfs,&pf,NULL,&kmax,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		Reduceloadx(pf,Kfs,ys);   xdelete(&Kfs);
Index: /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp	(revision 13876)
+++ /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp	(revision 13877)
@@ -24,5 +24,5 @@
 	Vector<IssmDouble>* ys  = NULL;
 
-	Loads* loads=NULL;
+	Loads* savedloads=NULL;
 	bool converged;
 	int constraints_converged;
@@ -33,5 +33,5 @@
 	int min_mechanical_constraints;
 	int max_nonlinear_iterations;
-	int  configuration_type;
+	int configuration_type;
 
 	/*Recover parameters: */
@@ -42,6 +42,7 @@
 
 	/*Were loads requested as output? : */
-	if(conserve_loads) loads=static_cast<Loads*>(femmodel->loads->Copy()); //protect loads from being modified by the solution
-	else               loads=static_cast<Loads*>(femmodel->loads);         //modify loads  in this solution
+	if(conserve_loads){
+		savedloads = static_cast<Loads*>(femmodel->loads->Copy());
+	}
 
 	count=1;
@@ -49,5 +50,5 @@
 
 	/*Start non-linear iteration using input velocity: */
-	GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices, loads, femmodel->materials, femmodel->parameters);
+	GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
 	Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
 
@@ -62,5 +63,5 @@
 		xdelete(&old_uf);old_uf=uf;
 
-		SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
+		femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		Reduceloadx(pf, Kfs, ys); xdelete(&Kfs);
@@ -72,5 +73,5 @@
 		InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
 
-		ConstraintsStatex(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
+		ConstraintsStatex(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 		if(VerboseConvergence()) _pprintLine_("   number of unstable constraints: " << num_unstable_constraints);
 
@@ -108,5 +109,8 @@
 
 	/*clean-up*/
-	if(conserve_loads) delete loads;
+	if(conserve_loads){
+		delete femmodel->loads;
+		femmodel->loads=savedloads;
+	}
 	xdelete(&uf);
 	xdelete(&ug);
Index: /issm/trunk-jpl/src/c/solvers/solver_stokescoupling_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_stokescoupling_nonlinear.cpp	(revision 13876)
+++ /issm/trunk-jpl/src/c/solvers/solver_stokescoupling_nonlinear.cpp	(revision 13877)
@@ -63,5 +63,5 @@
 
 		/*solve: */
-		SystemMatricesx(&Kff_horiz, &Kfs_horiz, &pf_horiz, &df_horiz, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		femmodel->SystemMatricesx(&Kff_horiz, &Kfs_horiz, &pf_horiz, &df_horiz, NULL);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		Reduceloadx(pf_horiz, Kfs_horiz, ys); xdelete(&Kfs_horiz);
@@ -77,5 +77,5 @@
 
 		/*solve: */
-		SystemMatricesx(&Kff_vert, &Kfs_vert, &pf_vert,  &df_vert,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		femmodel->SystemMatricesx(&Kff_vert, &Kfs_vert, &pf_vert,  &df_vert,NULL);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		Reduceloadx(pf_vert, Kfs_vert, ys); xdelete(&Kfs_vert);
Index: /issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp	(revision 13876)
+++ /issm/trunk-jpl/src/c/solvers/solver_thermal_nonlinear.cpp	(revision 13877)
@@ -53,5 +53,5 @@
 	for(;;){
 
-		SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		Reduceloadx(pf, Kfs, ys); xdelete(&Kfs); xdelete(&tf);