Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 13877)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 13878)
@@ -337,4 +337,46 @@
 
 /*Modules:*/
+void FemModel::AllocateSystemMatrices(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,Vector<IssmDouble>** pdf,Vector<IssmDouble>** ppf){ /*{{{*/
+
+	/*Intermediary*/
+	int      fsize,ssize;
+	int      connectivity, numberofdofspernode;
+	int      analysis_type, configuration_type;
+
+	/*output*/
+	Matrix<IssmDouble> *Kff  = NULL;
+	Matrix<IssmDouble> *Kfs  = NULL;
+	Vector<IssmDouble> *pf   = NULL;
+	Vector<IssmDouble> *df   = NULL;
+
+	bool oldalloc=true;
+
+	/*retrieve parameters: */
+	this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+	this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	this->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
+
+	/*retrieve node info*/
+	fsize=this->nodes->NumberOfDofs(configuration_type,FsetEnum);
+	ssize=this->nodes->NumberOfDofs(configuration_type,SsetEnum);
+	numberofdofspernode=this->nodes->MaxNumDofs(configuration_type,GsetEnum);
+
+	if(oldalloc){
+		Kff=new Matrix<IssmDouble>(fsize,fsize,connectivity,numberofdofspernode);
+		Kfs=new Matrix<IssmDouble>(fsize,ssize,connectivity,numberofdofspernode);
+		df =new Vector<IssmDouble>(fsize);
+		pf =new Vector<IssmDouble>(fsize);
+	}
+	else{
+		/*IN PROGRESS*/
+	}
+
+	/*Allocate output pointers*/
+	*pKff = Kff;
+	*pKfs = Kfs;
+	*pdf  = df;
+	*ppf  = pf;
+}
+/*}}}*/
 int  FemModel::UpdateVertexPositionsx(void){ /*{{{*/
 
@@ -522,7 +564,5 @@
 	/*intermediary: */
 	int      i;
-	int      fsize,ssize;
-	int      connectivity, numberofdofspernode;
-	int      analysis_type, configuration_type;
+	int      configuration_type;
 	Element *element = NULL;
 	Load    *load    = NULL;
@@ -535,30 +575,19 @@
 	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: */
+	this->AllocateSystemMatrices(&Kff,&Kfs,&df,&pf);
+
+	/*Fill stiffness matrix from elements and loads */
 	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));
@@ -566,40 +595,28 @@
 	}
 
-	/*Assemble matrices*/
+
+	/*Fill right hand side vector, from elements and loads */
+	for (i=0;i<this->elements->Size();i++){
+		element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		element->CreatePVector(pf);
+	}
+	for (i=0;i<this->loads->Size();i++){
+		load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
+		if (load->InAnalysis(configuration_type)) load->CreatePVector(pf);
+	}
+
+	/*Assemble matrices (required to get kmax)*/
 	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: */
+	/*Now, figure out maximum value of stiffness matrix, so that we can penalize it correctly: */
 	kmax=Kff->Norm(NORM_INF);
 
-	/*Now, deal with penalties*/
-	/*Fill stiffness matrix from loads: */
+	/*Now that we have kmax, deal with penalties (only in 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));
@@ -607,5 +624,7 @@
 	}
 
-	/*Assemble vector*/
+	/*Assemble matrices and vector*/
+	Kff->Assemble();
+	Kfs->Assemble();
 	pf->Assemble();
 
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 13877)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 13878)
@@ -52,4 +52,5 @@
 
 		/*Methods:*/
+		void AllocateSystemMatrices(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,Vector<IssmDouble>** pdf,Vector<IssmDouble>** ppf);
 		void Echo();
 		void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,const int* analyses,const int nummodels);