Index: /issm/trunk/src/c/parallel/thermal.cpp
===================================================================
--- /issm/trunk/src/c/parallel/thermal.cpp	(revision 653)
+++ /issm/trunk/src/c/parallel/thermal.cpp	(revision 654)
@@ -89,5 +89,4 @@
 	param=(Param*)femmodels[1].parameters->FindParamObject("p_g");
 	femmodels[1].parameters->DeleteObject((Object*)param);
-
 	
 	_printf_("call computational core:\n");
Index: /issm/trunk/src/c/parallel/thermal_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/thermal_core.cpp	(revision 653)
+++ /issm/trunk/src/c/parallel/thermal_core.cpp	(revision 654)
@@ -26,13 +26,13 @@
 
 	/*solutions vectors: */
-	Vec* t_g=NULL;
-	Vec* m_g=NULL;
+	Vec*    t_g=NULL;
+	Vec*    m_g=NULL;
 	double* time=NULL;
 
 	/*flags: */
-	int debug=0;
-	int numberofdofspernode;
-	int numberofnodes;
-	int nsteps;
+	int     debug=0;
+	int     numberofdofspernode;
+	int     numberofnodes;
+	int     nsteps;
 	double  dt;
 	double  ndt;
Index: /issm/trunk/src/c/parallel/thermal_core_nonlinear.cpp
===================================================================
--- /issm/trunk/src/c/parallel/thermal_core_nonlinear.cpp	(revision 653)
+++ /issm/trunk/src/c/parallel/thermal_core_nonlinear.cpp	(revision 654)
@@ -20,7 +20,9 @@
 	/*intermediary: */
 	Mat Kgg=NULL;
+	Mat Kgg_nopenalty=NULL;
 	Mat Kff=NULL;
 	Mat Kfs=NULL;
 	Vec pg=NULL;
+	Vec pg_nopenalty=NULL;
 	Vec pf=NULL;
 
@@ -36,4 +38,5 @@
 	char* solver_string=NULL;
 	int debug=0;
+	int lowmem=0;
 
 	/*Recover parameters: */
@@ -45,4 +48,5 @@
 	fem->parameters->FindParam((void*)&solver_string,"solverstring");
 	fem->parameters->FindParam((void*)&debug,"debug");
+	fem->parameters->FindParam((void*)&lowmem,"lowmem");
 	fem->parameters->FindParam((void*)&min_thermal_constraints,"min_thermal_constraints");
 
@@ -57,7 +61,25 @@
 
 		//*Generate system matrices
-		SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); 
-		//apply penalties
-		PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,inputs,analysis_type,sub_analysis_type); 
+		if (!lowmem){
+
+			/*Compute Kgg_nopenalty and pg_nopenalty once for all: */
+			if (count==1){
+				SystemMatricesx(&Kgg_nopenalty, &pg_nopenalty,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); 
+			}
+
+			/*Copy K_gg_nopenalty into Kgg, same for pg: */
+			Kgg=(Mat)xmalloc(sizeof(Mat));
+			MatDuplicate(Kgg_nopenalty,MAT_COPY_VALUES,&Kgg);
+			pg=(Vec)xmalloc(sizeof(Vec));
+			VecDuplicate(pg_nopenalty,&pg);VecCopy(pg_nopenalty,pg);
+
+			//apply penalties each time
+			PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,inputs,analysis_type,sub_analysis_type); 
+		}
+		else{
+			SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type); 
+			//apply penalties
+			PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,inputs,analysis_type,sub_analysis_type); 
+		}
 
 		/*!Reduce matrix from g to f size:*/