Index: /issm/trunk/src/c/parallel/thermal_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/thermal_core.cpp	(revision 617)
+++ /issm/trunk/src/c/parallel/thermal_core.cpp	(revision 618)
@@ -19,7 +19,9 @@
 
 	/*intermediary: */
+	Mat Kgg_nopenalty=NULL;
 	Mat Kgg=NULL;
 	Mat Kff=NULL;
 	Mat Kfs=NULL;
+	Vec pg_nopenalty=NULL;
 	Vec pg=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:*/