Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 16125)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 16126)
@@ -279,4 +279,10 @@
 					./modules/SurfaceAreax/SurfaceAreax.h\
 					./modules/SurfaceAreax/SurfaceAreax.cpp\
+					./modules/AllocateSystemMatricesx/AllocateSystemMatricesx.h\
+					./modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp\
+					./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h\
+					./modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp\
+					./modules/SystemMatricesx/SystemMatricesx.h\
+					./modules/SystemMatricesx/SystemMatricesx.cpp\
 					./modules/CreateNodalConstraintsx/CreateNodalConstraintsx.h\
 					./modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp\
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 16125)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 16126)
@@ -366,298 +366,4 @@
 /*}}}*/
 /*Modules:*/
-void FemModel::AllocateSystemMatrices(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,Vector<IssmDouble>** pdf,Vector<IssmDouble>** ppf){ /*{{{*/
-
-	/*Intermediary*/
-	int  fsize,ssize,flocalsize,slocalsize;
-	int  connectivity, numberofdofspernode;
-	int  configuration_type;
-	int  m,n,M,N;
-	int *d_nnz = NULL;
-	int *o_nnz = NULL;
-
-	/*output*/
-	Matrix<IssmDouble> *Kff  = NULL;
-	Matrix<IssmDouble> *Kfs  = NULL;
-	Vector<IssmDouble> *pf   = NULL;
-	Vector<IssmDouble> *df   = NULL;
-
-	bool oldalloc=false;
-
-	/*retrieve parameters: */
-	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);
-	flocalsize = this->nodes->NumberOfDofsLocal(configuration_type,FsetEnum);
-	slocalsize = this->nodes->NumberOfDofsLocal(configuration_type,SsetEnum);
-
-	numberofdofspernode=this->nodes->MaxNumDofs(configuration_type,GsetEnum);
-
-	if(oldalloc){
-		if(pKff) Kff=new Matrix<IssmDouble>(fsize,fsize,connectivity,numberofdofspernode);
-		if(pKfs) Kfs=new Matrix<IssmDouble>(fsize,ssize,connectivity,numberofdofspernode);
-		if(pdf)  df =new Vector<IssmDouble>(fsize);
-		if(ppf)  pf =new Vector<IssmDouble>(fsize);
-	}
-	else{
-		if(pKff){
-			m=flocalsize; n=flocalsize; /*local  sizes*/
-			M=fsize;      N=fsize;      /*global sizes*/
-			this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,FsetEnum);
-			Kff=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
-			xDelete<int>(d_nnz);
-			xDelete<int>(o_nnz);
-		}
-		if(pKfs){
-			m=flocalsize; n=slocalsize; /*local  sizes*/
-			M=fsize;      N=ssize;      /*global sizes*/
-			this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,SsetEnum);
-			Kfs=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
-			xDelete<int>(d_nnz);
-			xDelete<int>(o_nnz);
-		}
-		if(pdf) df =new Vector<IssmDouble>(flocalsize,fsize);
-		if(ppf) pf =new Vector<IssmDouble>(flocalsize,fsize);
-	}
-
-	/*Allocate output pointers*/
-	if(pKff) *pKff = Kff;
-	if(pKfs) *pKfs = Kfs;
-	if(pdf)  *pdf  = df;
-	if(ppf)  *ppf  = pf;
-}
-/*}}}*/
-void FemModel::MatrixNonzeros(int** pd_nnz,int** po_nnz,int set1enum,int set2enum){/*{{{*/
-
-	/*Intermediary*/
-	int      i,j,k,index,offset,count;
-	int      configuration_type;
-	int      d_nz,o_nz;
-	Element *element            = NULL;
-	Load    *load               = NULL;
-	int     *head_e             = NULL;
-	int     *next_e             = NULL;
-	int     *count2offset_e     = NULL;
-	int     *head_l             = NULL;
-	int     *next_l             = NULL;
-	int     *count2offset_l     = NULL;
-	int     *lidlist            = NULL;
-
-	/*output*/
-	int *d_nnz = NULL;
-	int *o_nnz = NULL;
-
-	/*retrive parameters: */
-	this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
-
-	/*Get vector size and number of nodes*/
-	int numnodes            = nodes->NumberOfNodes(configuration_type);
-	int localnumnodes       = nodes->Size();
-	int numberofdofspernode = nodes->MaxNumDofs(configuration_type,GsetEnum);
-	int M                   = nodes->NumberOfDofs(configuration_type,set1enum);
-	int N                   = nodes->NumberOfDofs(configuration_type,set2enum);
-	int m                   = nodes->NumberOfDofsLocal(configuration_type,set1enum);
-	int n                   = nodes->NumberOfDofsLocal(configuration_type,set2enum);
-	int numnodesperelement  = elements->MaxNumNodes();
-	int numnodesperload     = loads->MaxNumNodes(configuration_type);
-
-	/*First, we are building chaining vectors so that we know what nodes are
-	 * connected to what elements. These vectors are such that:
-	 *   for(int i=head[id];i!=-1;i=next[i])
-	 * will loop over all the elements that are connected to the node number
-	 * id*/
-	head_e         = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_e[i]=-1;
-	next_e         = xNew<int>(elements->Size()*numnodesperelement);
-	count2offset_e = xNew<int>(elements->Size()*numnodesperelement);
-
-	k=0;
-	for(i=0;i<elements->Size();i++){
-		element = dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-		lidlist = xNew<int>(element->GetNumberOfNodes());
-		element->GetNodesLidList(lidlist);
-
-		for(j=0;j<element->GetNumberOfNodes();j++){
-			index = lidlist[j];
-			_assert_(index>=0 && index<numnodes);
-
-			count2offset_e[k]=i;
-			next_e[k]=head_e[index];
-			head_e[index]=k++;
-		}
-		for(j=0;j<numnodesperelement-element->GetNumberOfNodes();j++) k++;
-
-		xDelete<int>(lidlist);
-	}
-
-	/*Chain for loads*/
-	head_l         = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_l[i]=-1;
-	next_l         = xNew<int>(loads->Size(configuration_type)*numnodesperload);
-	count2offset_l = xNew<int>(loads->Size(configuration_type)*numnodesperload);
-	k=0;
-	for(i=0;i<loads->Size();i++){
-		load = dynamic_cast<Load*>(loads->GetObjectByOffset(i));
-		if(!load->InAnalysis(configuration_type)) continue;
-		lidlist = xNew<int>(load->GetNumberOfNodes());
-		load->GetNodesLidList(lidlist);
-
-		for(j=0;j<load->GetNumberOfNodes();j++){
-			index = lidlist[j];
-			_assert_(index>=0 && index<numnodes);
-
-			count2offset_l[k]=i;
-			next_l[k]=head_l[index];
-			head_l[index]=k++;
-		}
-		for(j=0;j<numnodesperload-load->GetNumberOfNodes();j++) k++;
-
-		xDelete<int>(lidlist);
-	}
-
-	/*OK now count number of dofs and flag each nodes for each node i*/
-	bool *flags                  = xNew<bool>(localnumnodes);
-	int  *flagsindices           = xNew<int>(localnumnodes);
-	int  *d_connectivity         = xNewZeroInit<int>(numnodes);
-	int  *o_connectivity         = xNewZeroInit<int>(numnodes);
-	int  *connectivity_clone     = xNewZeroInit<int>(numnodes);
-	int  *all_connectivity_clone = xNewZeroInit<int>(numnodes);
-
-	/*Resetting flags to false at eahc iteration takes a lot of time, so we keep track of the flags
-	 * to reset in flagsindices, initialized with -1*/
-	for(i = 0;i<localnumnodes;i++) flags[i]        = false;
-	for(i = 0;i<localnumnodes;i++) flagsindices[i] = -1;
-
-	/*Create connectivity vector*/
-	for(i=0;i<nodes->Size();i++){
-		Node* node=dynamic_cast<Node*>(nodes->GetObjectByOffset(i));
-		if(node->InAnalysis(configuration_type)){
-
-			/*Reinitialize flags to false*/
-			j=0;
-			while(true){
-				if(flagsindices[j]>=0){
-					flags[flagsindices[j]] = false;
-					flagsindices[j]        = -1;
-					j++;
-				}
-				else{
-					break;
-				}
-			}
-
-			//for(j=0;j<localnumnodes;j++) flags[j]=false;
-
-			/*Loop over elements that hold node number i*/
-			//if(head_e[node->Lid()]==-1 && head_l[node->Lid()]==-1){
-			//	printf("[%i] vertex %i\n",IssmComm::GetRank(),node->Lid()+1);
-			//}
-			for(j=head_e[node->Lid()];j!=-1;j=next_e[j]){
-				offset=count2offset_e[j];
-				element=dynamic_cast<Element*>(elements->GetObjectByOffset(offset));
-				element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
-				if(node->IsClone()){
-					connectivity_clone[node->Sid()]+=d_nz+o_nz;
-				}
-				else{
-					d_connectivity[node->Sid()]+=d_nz;
-					o_connectivity[node->Sid()]+=o_nz;
-				}
-			}
-			for(j=head_l[node->Lid()];j!=-1;j=next_l[j]){
-				offset=count2offset_l[j];
-				load=dynamic_cast<Load*>(loads->GetObjectByOffset(offset));
-				load->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
-				if(node->IsClone()){
-					connectivity_clone[node->Sid()]+=d_nz+o_nz;
-				}
-				else{
-					d_connectivity[node->Sid()]+=d_nz;
-					o_connectivity[node->Sid()]+=o_nz;
-				}
-			}
-		}
-	}
-	xDelete<bool>(flags);
-	xDelete<int>(flagsindices);
-	xDelete<int>(count2offset_e);
-	xDelete<int>(head_e);
-	xDelete<int>(next_e);
-	xDelete<int>(count2offset_l);
-	xDelete<int>(head_l);
-	xDelete<int>(next_l);
-
-	/*sum over all cpus*/
-	ISSM_MPI_Allreduce((void*)connectivity_clone,(void*)all_connectivity_clone,numnodes,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm());
-	xDelete<int>(connectivity_clone);
-
-	if(set1enum==FsetEnum){
-		count=0;
-		d_nnz=xNew<int>(m);
-		o_nnz=xNew<int>(m);
-		for(i=0;i<nodes->Size();i++){
-			Node* node=dynamic_cast<Node*>(nodes->GetObjectByOffset(i));
-			if(node->InAnalysis(configuration_type) && !node->IsClone()){
-				for(j=0;j<node->indexing.fsize;j++){
-					_assert_(count<m);
-					d_nnz[count]=numberofdofspernode*(d_connectivity[node->Sid()] + all_connectivity_clone[node->Sid()]);
-					o_nnz[count]=numberofdofspernode*(o_connectivity[node->Sid()] + all_connectivity_clone[node->Sid()]);
-					if(d_nnz[count]>n)   d_nnz[count]=n;
-					if(o_nnz[count]>N-n) o_nnz[count]=N-n;
-					count++;
-				}
-			}
-		}
-		_assert_(m==count);
-	}
-	else{
-		_error_("STOP not implemented");
-	}
-	xDelete<int>(d_connectivity);
-	xDelete<int>(o_connectivity);
-	xDelete<int>(all_connectivity_clone);
-
-	/*Allocate ouptput pointer*/
-	*pd_nnz=d_nnz;
-	*po_nnz=o_nnz;
-
-}/*}}}*/
-void FemModel::CreateJacobianMatrixx(Matrix<IssmDouble>** pJff,IssmDouble kmax){/*{{{*/
-
-	int      i,connectivity;
-	int      numberofdofspernode;
-	int      fsize,configuration_type;
-	Element *element = NULL;
-	Load    *load    = NULL;
-	Matrix<IssmDouble>* Jff = NULL;
-
-	/*Checks*/
-	_assert_(nodes && elements);
-
-	/*Recover some parameters*/
-	parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
-	parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
-	fsize=nodes->NumberOfDofs(configuration_type,FsetEnum);
-	numberofdofspernode=nodes->MaxNumDofs(configuration_type,GsetEnum);
-
-	/*Initialize Jacobian Matrix*/
-	this->AllocateSystemMatrices(&Jff,NULL,NULL,NULL);
-
-	/*Create and assemble matrix*/
-	for(i=0;i<elements->Size();i++){
-		element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-		element->CreateJacobianMatrix(Jff);
-	}
-	for (i=0;i<loads->Size();i++){
-		load=(Load*)loads->GetObjectByOffset(i);
-		if(load->InAnalysis(configuration_type)) load->CreateJacobianMatrix(Jff);
-		if(load->InAnalysis(configuration_type)) load->PenaltyCreateJacobianMatrix(Jff,kmax);
-	}
-	Jff->Assemble();
-
-	/*Assign output pointer*/
-	*pJff=Jff;
-
-}/*}}}*/
 int  FemModel::UpdateVertexPositionsx(void){ /*{{{*/
 
@@ -845,117 +551,4 @@
 }
 /*}}}*/
-void FemModel::SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax){/*{{{*/
-
-	/*intermediary: */
-	int      i,M,N;
-	int      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 = 0;
-
-	/*Display message*/
-	if(VerboseModule()) _printf0_("   Generating matrices");
-
-	/*retrive parameters: */
-	this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
-
-	/*First, we might need to do a dry run to get kmax if penalties are employed*/
-	if(loads->IsPenalty(configuration_type)){
-
-		/*Allocate Kff_temp*/
-		Matrix<IssmDouble> *Kff_temp = NULL;
-		this->AllocateSystemMatrices(&Kff_temp,NULL,NULL,NULL);
-
-		/*Get complete stiffness matrix without penalties*/
-		for (i=0;i<this->elements->Size();i++){
-			element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
-			element->CreateKMatrix(Kff_temp,NULL);
-		}
-
-		for (i=0;i<this->loads->Size();i++){
-			load=dynamic_cast<Load*>(this->loads->GetObjectByOffset(i));
-			if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff_temp,NULL);
-		}
-		Kff_temp->Assemble();
-
-		/*Now, figure out maximum value of stiffness matrix, so that we can penalize it correctly: */
-		kmax=Kff_temp->Norm(NORM_INF);
-		delete Kff_temp;
-	}
-
-	/*Allocate stiffness matrices and load vector*/
-	this->AllocateSystemMatrices(&Kff,&Kfs,&df,&pf);
-
-	/*Display size*/
-	if(VerboseModule()){
-		Kff->GetSize(&M,&N);
-		_printf0_(" (Kff stiffness matrix size: "<<M<<" x "<<N<<")\n");
-	}
-
-	/*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);
-	}
-
-	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);
-	}
-
-	/*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);
-	}
-
-	/*Now deal with penalties (only in loads)*/
-	if(loads->IsPenalty(configuration_type)){
-		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);
-		}
-		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);
-		}
-	}
-
-	/*Create dof vector for stiffness matrix preconditioning*/
-	for (i=0;i<this->elements->Size();i++){
-		element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
-		element->CreateDVector(df);
-	}
-
-	/*Assemble matrices and vector*/
-	Kff->Assemble();
-	Kfs->Assemble();
-	pf->Assemble();
-	df->Assemble();
-	//Kff->AllocationInfo();
-	//Kfs->AllocationInfo();
-
-	/*Assign output pointers: */
-	if(pKff) *pKff=Kff;
-	else      delete Kff;
-	if(pKfs) *pKfs=Kfs;
-	else      delete Kfs;
-	if(ppf)  *ppf=pf;
-	else      delete pf;
-	if(pdf)  *pdf=df;
-	else      delete 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 16125)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 16126)
@@ -51,9 +51,6 @@
 
 		/*Methods:*/
-		void AllocateSystemMatrices(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,Vector<IssmDouble>** pdf,Vector<IssmDouble>** ppf);
-		void CreateJacobianMatrixx(Matrix<IssmDouble>** pJff,IssmDouble kmax);
 		void Echo();
 		void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,const int* analyses,const int nummodels);
-		void MatrixNonzeros(int** pd_nnz,int** po_nnz,int set1enum,int set2enum);
 		void Solve(void);
 		void OutputResults(void);
@@ -95,5 +92,4 @@
 		void Deflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y);
 		#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/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp	(revision 16126)
+++ /issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp	(revision 16126)
@@ -0,0 +1,264 @@
+/*!\file AllocateSystemMatricesx
+ * \brief retrieve vector from inputs in elements
+ */
+
+#include "./AllocateSystemMatricesx.h"
+#include "../../shared/shared.h"
+#include "../../toolkits/toolkits.h"
+
+void AllocateSystemMatricesx(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,Vector<IssmDouble>** pdf,Vector<IssmDouble>** ppf,FemModel* femmodel){
+
+	/*Intermediary*/
+	int  fsize,ssize,flocalsize,slocalsize;
+	int  connectivity, numberofdofspernode;
+	int  configuration_type;
+	int  m,n,M,N;
+	int *d_nnz = NULL;
+	int *o_nnz = NULL;
+
+	/*output*/
+	Matrix<IssmDouble> *Kff  = NULL;
+	Matrix<IssmDouble> *Kfs  = NULL;
+	Vector<IssmDouble> *pf   = NULL;
+	Vector<IssmDouble> *df   = NULL;
+
+	bool oldalloc=false;
+
+	/*retrieve parameters: */
+	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	femmodel->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
+
+	/*retrieve node info*/
+	fsize      = femmodel->nodes->NumberOfDofs(configuration_type,FsetEnum);
+	ssize      = femmodel->nodes->NumberOfDofs(configuration_type,SsetEnum);
+	flocalsize = femmodel->nodes->NumberOfDofsLocal(configuration_type,FsetEnum);
+	slocalsize = femmodel->nodes->NumberOfDofsLocal(configuration_type,SsetEnum);
+
+	numberofdofspernode=femmodel->nodes->MaxNumDofs(configuration_type,GsetEnum);
+
+	if(oldalloc){
+		if(pKff) Kff=new Matrix<IssmDouble>(fsize,fsize,connectivity,numberofdofspernode);
+		if(pKfs) Kfs=new Matrix<IssmDouble>(fsize,ssize,connectivity,numberofdofspernode);
+		if(pdf)  df =new Vector<IssmDouble>(fsize);
+		if(ppf)  pf =new Vector<IssmDouble>(fsize);
+	}
+	else{
+		if(pKff){
+			m=flocalsize; n=flocalsize; /*local  sizes*/
+			M=fsize;      N=fsize;      /*global sizes*/
+			MatrixNonzeros(&d_nnz,&o_nnz,femmodel,FsetEnum,FsetEnum);
+			Kff=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
+			xDelete<int>(d_nnz);
+			xDelete<int>(o_nnz);
+		}
+		if(pKfs){
+			m=flocalsize; n=slocalsize; /*local  sizes*/
+			M=fsize;      N=ssize;      /*global sizes*/
+			MatrixNonzeros(&d_nnz,&o_nnz,femmodel,FsetEnum,SsetEnum);
+			Kfs=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
+			xDelete<int>(d_nnz);
+			xDelete<int>(o_nnz);
+		}
+		if(pdf) df =new Vector<IssmDouble>(flocalsize,fsize);
+		if(ppf) pf =new Vector<IssmDouble>(flocalsize,fsize);
+	}
+
+	/*Allocate output pointers*/
+	if(pKff) *pKff = Kff;
+	if(pKfs) *pKfs = Kfs;
+	if(pdf)  *pdf  = df;
+	if(ppf)  *ppf  = pf;
+}
+
+void MatrixNonzeros(int** pd_nnz,int** po_nnz,FemModel* femmodel,int set1enum,int set2enum){
+
+	/*Intermediary*/
+	int      i,j,k,index,offset,count;
+	int      configuration_type;
+	int      d_nz,o_nz;
+	Element *element            = NULL;
+	Load    *load               = NULL;
+	int     *head_e             = NULL;
+	int     *next_e             = NULL;
+	int     *count2offset_e     = NULL;
+	int     *head_l             = NULL;
+	int     *next_l             = NULL;
+	int     *count2offset_l     = NULL;
+	int     *lidlist            = NULL;
+
+	/*output*/
+	int *d_nnz = NULL;
+	int *o_nnz = NULL;
+
+	/*retrive parameters: */
+	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+
+	/*Get vector size and number of nodes*/
+	int numnodes            = femmodel->nodes->NumberOfNodes(configuration_type);
+	int localnumnodes       = femmodel->nodes->Size();
+	int numberofdofspernode = femmodel->nodes->MaxNumDofs(configuration_type,GsetEnum);
+	int M                   = femmodel->nodes->NumberOfDofs(configuration_type,set1enum);
+	int N                   = femmodel->nodes->NumberOfDofs(configuration_type,set2enum);
+	int m                   = femmodel->nodes->NumberOfDofsLocal(configuration_type,set1enum);
+	int n                   = femmodel->nodes->NumberOfDofsLocal(configuration_type,set2enum);
+	int numnodesperelement  = femmodel->elements->MaxNumNodes();
+	int numnodesperload     = femmodel->loads->MaxNumNodes(configuration_type);
+
+	/*First, we are building chaining vectors so that we know what nodes are
+	 * connected to what elements. These vectors are such that:
+	 *   for(int i=head[id];i!=-1;i=next[i])
+	 * will loop over all the elements that are connected to the node number
+	 * id*/
+	head_e         = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_e[i]=-1;
+	next_e         = xNew<int>(femmodel->elements->Size()*numnodesperelement);
+	count2offset_e = xNew<int>(femmodel->elements->Size()*numnodesperelement);
+
+	k=0;
+	for(i=0;i<femmodel->elements->Size();i++){
+		element = dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		lidlist = xNew<int>(element->GetNumberOfNodes());
+		element->GetNodesLidList(lidlist);
+
+		for(j=0;j<element->GetNumberOfNodes();j++){
+			index = lidlist[j];
+			_assert_(index>=0 && index<numnodes);
+
+			count2offset_e[k]=i;
+			next_e[k]=head_e[index];
+			head_e[index]=k++;
+		}
+		for(j=0;j<numnodesperelement-element->GetNumberOfNodes();j++) k++;
+
+		xDelete<int>(lidlist);
+	}
+
+	/*Chain for loads*/
+	head_l         = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_l[i]=-1;
+	next_l         = xNew<int>(femmodel->loads->Size(configuration_type)*numnodesperload);
+	count2offset_l = xNew<int>(femmodel->loads->Size(configuration_type)*numnodesperload);
+	k=0;
+	for(i=0;i<femmodel->loads->Size();i++){
+		load = dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
+		if(!load->InAnalysis(configuration_type)) continue;
+		lidlist = xNew<int>(load->GetNumberOfNodes());
+		load->GetNodesLidList(lidlist);
+
+		for(j=0;j<load->GetNumberOfNodes();j++){
+			index = lidlist[j];
+			_assert_(index>=0 && index<numnodes);
+
+			count2offset_l[k]=i;
+			next_l[k]=head_l[index];
+			head_l[index]=k++;
+		}
+		for(j=0;j<numnodesperload-load->GetNumberOfNodes();j++) k++;
+
+		xDelete<int>(lidlist);
+	}
+
+	/*OK now count number of dofs and flag each nodes for each node i*/
+	bool *flags                  = xNew<bool>(localnumnodes);
+	int  *flagsindices           = xNew<int>(localnumnodes);
+	int  *d_connectivity         = xNewZeroInit<int>(numnodes);
+	int  *o_connectivity         = xNewZeroInit<int>(numnodes);
+	int  *connectivity_clone     = xNewZeroInit<int>(numnodes);
+	int  *all_connectivity_clone = xNewZeroInit<int>(numnodes);
+
+	/*Resetting flags to false at eahc iteration takes a lot of time, so we keep track of the flags
+	 * to reset in flagsindices, initialized with -1*/
+	for(i = 0;i<localnumnodes;i++) flags[i]        = false;
+	for(i = 0;i<localnumnodes;i++) flagsindices[i] = -1;
+
+	/*Create connectivity vector*/
+	for(i=0;i<femmodel->nodes->Size();i++){
+		Node* node=dynamic_cast<Node*>(femmodel->nodes->GetObjectByOffset(i));
+		if(node->InAnalysis(configuration_type)){
+
+			/*Reinitialize flags to false*/
+			j=0;
+			while(true){
+				if(flagsindices[j]>=0){
+					flags[flagsindices[j]] = false;
+					flagsindices[j]        = -1;
+					j++;
+				}
+				else{
+					break;
+				}
+			}
+
+			//for(j=0;j<localnumnodes;j++) flags[j]=false;
+
+			/*Loop over elements that hold node number i*/
+			//if(head_e[node->Lid()]==-1 && head_l[node->Lid()]==-1){
+			//	printf("[%i] vertex %i\n",IssmComm::GetRank(),node->Lid()+1);
+			//}
+			for(j=head_e[node->Lid()];j!=-1;j=next_e[j]){
+				offset=count2offset_e[j];
+				element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(offset));
+				element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
+				if(node->IsClone()){
+					connectivity_clone[node->Sid()]+=d_nz+o_nz;
+				}
+				else{
+					d_connectivity[node->Sid()]+=d_nz;
+					o_connectivity[node->Sid()]+=o_nz;
+				}
+			}
+			for(j=head_l[node->Lid()];j!=-1;j=next_l[j]){
+				offset=count2offset_l[j];
+				load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(offset));
+				load->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
+				if(node->IsClone()){
+					connectivity_clone[node->Sid()]+=d_nz+o_nz;
+				}
+				else{
+					d_connectivity[node->Sid()]+=d_nz;
+					o_connectivity[node->Sid()]+=o_nz;
+				}
+			}
+		}
+	}
+	xDelete<bool>(flags);
+	xDelete<int>(flagsindices);
+	xDelete<int>(count2offset_e);
+	xDelete<int>(head_e);
+	xDelete<int>(next_e);
+	xDelete<int>(count2offset_l);
+	xDelete<int>(head_l);
+	xDelete<int>(next_l);
+
+	/*sum over all cpus*/
+	ISSM_MPI_Allreduce((void*)connectivity_clone,(void*)all_connectivity_clone,numnodes,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm());
+	xDelete<int>(connectivity_clone);
+
+	if(set1enum==FsetEnum){
+		count=0;
+		d_nnz=xNew<int>(m);
+		o_nnz=xNew<int>(m);
+		for(i=0;i<femmodel->nodes->Size();i++){
+			Node* node=dynamic_cast<Node*>(femmodel->nodes->GetObjectByOffset(i));
+			if(node->InAnalysis(configuration_type) && !node->IsClone()){
+				for(j=0;j<node->indexing.fsize;j++){
+					_assert_(count<m);
+					d_nnz[count]=numberofdofspernode*(d_connectivity[node->Sid()] + all_connectivity_clone[node->Sid()]);
+					o_nnz[count]=numberofdofspernode*(o_connectivity[node->Sid()] + all_connectivity_clone[node->Sid()]);
+					if(d_nnz[count]>n)   d_nnz[count]=n;
+					if(o_nnz[count]>N-n) o_nnz[count]=N-n;
+					count++;
+				}
+			}
+		}
+		_assert_(m==count);
+	}
+	else{
+		_error_("STOP not implemented");
+	}
+	xDelete<int>(d_connectivity);
+	xDelete<int>(o_connectivity);
+	xDelete<int>(all_connectivity_clone);
+
+	/*Allocate ouptput pointer*/
+	*pd_nnz=d_nnz;
+	*po_nnz=o_nnz;
+}
Index: /issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.h	(revision 16126)
+++ /issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.h	(revision 16126)
@@ -0,0 +1,14 @@
+/*!\file:  AllocateSystemMatricesx.h
+*/ 
+
+#ifndef _ALLOCATESYSTEMMATRICESX_H
+#define _ALLOCATESYSTEMMATRICESX_H
+
+#include "../../classes/classes.h"
+
+/* local prototypes: */
+void AllocateSystemMatricesx(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,Vector<IssmDouble>** pdf,Vector<IssmDouble>** ppf,FemModel* femmodel);
+void MatrixNonzeros(int** pd_nnz,int** po_nnz,FemModel* femmodel,int set1enum,int set2enum);
+
+
+#endif
Index: /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp	(revision 16126)
+++ /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.cpp	(revision 16126)
@@ -0,0 +1,46 @@
+/*!\file CreateJacobianMatrixx
+ * \brief retrieve vector from inputs in elements
+ */
+
+#include "./CreateJacobianMatrixx.h"
+#include "../../shared/shared.h"
+#include "../../toolkits/toolkits.h"
+#include "../AllocateSystemMatricesx/AllocateSystemMatricesx.h"
+
+void CreateJacobianMatrixx(Matrix<IssmDouble>** pJff,FemModel* femmodel,IssmDouble kmax){
+
+	int      i,connectivity;
+	int      numberofdofspernode;
+	int      fsize,configuration_type;
+	Element *element = NULL;
+	Load    *load    = NULL;
+	Matrix<IssmDouble>* Jff = NULL;
+
+	/*Checks*/
+	_assert_(femmodel && femmodel->nodes && femmodel->elements);
+
+	/*Recover some parameters*/
+	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	femmodel->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
+	fsize=femmodel->nodes->NumberOfDofs(configuration_type,FsetEnum);
+	numberofdofspernode=femmodel->nodes->MaxNumDofs(configuration_type,GsetEnum);
+
+	/*Initialize Jacobian Matrix*/
+	AllocateSystemMatricesx(&Jff,NULL,NULL,NULL,femmodel);
+
+	/*Create and assemble matrix*/
+	for(i=0;i<femmodel->elements->Size();i++){
+		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		element->CreateJacobianMatrix(Jff);
+	}
+	for (i=0;i<femmodel->loads->Size();i++){
+		load=(Load*)femmodel->loads->GetObjectByOffset(i);
+		if(load->InAnalysis(configuration_type)) load->CreateJacobianMatrix(Jff);
+		if(load->InAnalysis(configuration_type)) load->PenaltyCreateJacobianMatrix(Jff,kmax);
+	}
+	Jff->Assemble();
+
+	/*Assign output pointer*/
+	*pJff=Jff;
+
+}
Index: /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h	(revision 16126)
+++ /issm/trunk-jpl/src/c/modules/CreateJacobianMatrixx/CreateJacobianMatrixx.h	(revision 16126)
@@ -0,0 +1,12 @@
+/*!\file:  CreateJacobianMatrixx.h
+*/ 
+
+#ifndef _CREATEJACOBIANMATRIXX_H
+#define _CREATEJACOBIANMATRIXX_H
+
+#include "../../classes/classes.h"
+
+/* local prototypes: */
+void CreateJacobianMatrixx(Matrix<IssmDouble>** pJff,FemModel* femmodel,IssmDouble kmax);
+
+#endif
Index: /issm/trunk-jpl/src/c/modules/Exp2Kmlx/Exp2Kmlx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/Exp2Kmlx/Exp2Kmlx.h	(revision 16125)
+++ /issm/trunk-jpl/src/c/modules/Exp2Kmlx/Exp2Kmlx.h	(revision 16126)
@@ -10,10 +10,6 @@
 
 /* local prototypes: */
-int Exp2Kmlx(char* filexp,char* filkml,
-			 int sgn,
-			 bool holes);
-int Exp2Kmlx(char* filexp,char* filkml,
-			 int sgn,double cm,double sp,
-			 bool holes);
+int Exp2Kmlx(char* filexp,char* filkml,int sgn,bool holes);
+int Exp2Kmlx(char* filexp,char* filkml,int sgn,double cm,double sp,bool holes);
 
 #endif  /* _EXP2KMLX_H */
Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 16126)
+++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 16126)
@@ -0,0 +1,121 @@
+/*!\file SystemMatricesx
+ * \brief retrieve vector from inputs in elements
+ */
+
+#include "./SystemMatricesx.h"
+#include "../../shared/shared.h"
+#include "../../toolkits/toolkits.h"
+#include "../AllocateSystemMatricesx/AllocateSystemMatricesx.h"
+
+void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax,FemModel* femmodel){
+
+	/*intermediary: */
+	int      i,M,N;
+	int      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 = 0;
+
+	/*Display message*/
+	if(VerboseModule()) _printf0_("   Generating matrices");
+
+	/*retrive parameters: */
+	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+
+	/*First, we might need to do a dry run to get kmax if penalties are employed*/
+	if(femmodel->loads->IsPenalty(configuration_type)){
+
+		/*Allocate Kff_temp*/
+		Matrix<IssmDouble> *Kff_temp = NULL;
+		AllocateSystemMatricesx(&Kff_temp,NULL,NULL,NULL,femmodel);
+
+		/*Get complete stiffness matrix without penalties*/
+		for (i=0;i<femmodel->elements->Size();i++){
+			element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+			element->CreateKMatrix(Kff_temp,NULL);
+		}
+
+		for (i=0;i<femmodel->loads->Size();i++){
+			load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
+			if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff_temp,NULL);
+		}
+		Kff_temp->Assemble();
+
+		/*Now, figure out maximum value of stiffness matrix, so that we can penalize it correctly: */
+		kmax=Kff_temp->Norm(NORM_INF);
+		delete Kff_temp;
+	}
+
+	/*Allocate stiffness matrices and load vector*/
+	AllocateSystemMatricesx(&Kff,&Kfs,&df,&pf,femmodel);
+
+	/*Display size*/
+	if(VerboseModule()){
+		Kff->GetSize(&M,&N);
+		_printf0_(" (Kff stiffness matrix size: "<<M<<" x "<<N<<")\n");
+	}
+
+	/*Fill stiffness matrix from elements and loads */
+	for (i=0;i<femmodel->elements->Size();i++){
+		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		element->CreateKMatrix(Kff,Kfs);
+	}
+
+	for (i=0;i<femmodel->loads->Size();i++){
+		load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
+		if(load->InAnalysis(configuration_type)) load->CreateKMatrix(Kff,Kfs);
+	}
+
+	/*Fill right hand side vector, from elements and loads */
+	for (i=0;i<femmodel->elements->Size();i++){
+		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		element->CreatePVector(pf);
+	}
+	for (i=0;i<femmodel->loads->Size();i++){
+		load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
+		if(load->InAnalysis(configuration_type)) load->CreatePVector(pf);
+	}
+
+	/*Now deal with penalties (only in loads)*/
+	if(femmodel->loads->IsPenalty(configuration_type)){
+		for (i=0;i<femmodel->loads->Size();i++){
+			load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
+			if(load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kff,Kfs,kmax);
+		}
+		for (i=0;i<femmodel->loads->Size();i++){
+			load=dynamic_cast<Load*>(femmodel->loads->GetObjectByOffset(i));
+			if(load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pf,kmax);
+		}
+	}
+
+	/*Create dof vector for stiffness matrix preconditioning*/
+	for (i=0;i<femmodel->elements->Size();i++){
+		element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		element->CreateDVector(df);
+	}
+
+	/*Assemble matrices and vector*/
+	Kff->Assemble();
+	Kfs->Assemble();
+	pf->Assemble();
+	df->Assemble();
+	//Kff->AllocationInfo();
+	//Kfs->AllocationInfo();
+
+	/*Assign output pointers: */
+	if(pKff) *pKff=Kff;
+	else      delete Kff;
+	if(pKfs) *pKfs=Kfs;
+	else      delete Kfs;
+	if(ppf)  *ppf=pf;
+	else      delete pf;
+	if(pdf)  *pdf=df;
+	else      delete df;
+	if(pkmax) *pkmax=kmax;
+}
Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h	(revision 16126)
+++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h	(revision 16126)
@@ -0,0 +1,12 @@
+/*!\file:  SystemMatricesx.h
+*/ 
+
+#ifndef _SYSTEMMATRICESX_H
+#define _SYSTEMMATRICESX_H
+
+#include "../../classes/classes.h"
+
+/* local prototypes: */
+void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax,FemModel* femmodel);
+
+#endif  /* _SYSTEMMATRICESX_H */
Index: /issm/trunk-jpl/src/c/modules/modules.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/modules.h	(revision 16125)
+++ /issm/trunk-jpl/src/c/modules/modules.h	(revision 16126)
@@ -7,4 +7,5 @@
 
 /*Modules: */
+#include "./AllocateSystemMatricesx/AllocateSystemMatricesx.h"
 #include "./AverageFilterx/AverageFilterx.h"
 #include "./AverageOntoPartitionx/AverageOntoPartitionx.h"
@@ -22,4 +23,5 @@
 #include "./ControlInputScaleGradientx/ControlInputScaleGradientx.h"
 #include "./CreateNodalConstraintsx/CreateNodalConstraintsx.h"
+#include "./CreateJacobianMatrixx/CreateJacobianMatrixx.h"
 #include "./Delta18oParameterizationx/Delta18oParameterizationx.h"
 #include "./DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h"
@@ -88,4 +90,5 @@
 #include "./SmbGradientsx/SmbGradientsx.h"
 #include "./Solverx/Solverx.h"
+#include "./SystemMatricesx/SystemMatricesx.h"
 #include "./SpcNodesx/SpcNodesx.h"
 #include "./SurfaceAreax/SurfaceAreax.h"
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_adjoint_linear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_adjoint_linear.cpp	(revision 16125)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_adjoint_linear.cpp	(revision 16126)
@@ -26,5 +26,5 @@
 	femmodel->UpdateConstraintsx();
 
-	femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
+	SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL,femmodel);
 	CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 	Reduceloadx(pf, Kfs, ys,true); delete Kfs; //true means spc = 0
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp	(revision 16125)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp	(revision 16126)
@@ -79,5 +79,5 @@
 		for(;;){
 			femmodel->HydrologyTransferx();
-			femmodel->SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax);
+			SystemMatricesx(&Kff,&Kfs,&pf,&df,&sediment_kmax,femmodel);
 			CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum);
 			Reduceloadx(pf,Kfs,ys); delete Kfs;
@@ -123,5 +123,5 @@
 			for(;;){
 				femmodel->HydrologyTransferx();
-				femmodel->SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL);
+				SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
 				CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);
 				Reduceloadx(pf,Kfs,ys); delete Kfs;
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp	(revision 16125)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp	(revision 16126)
@@ -24,5 +24,5 @@
 	femmodel->UpdateConstraintsx();
 
-	femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
+	SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
 	CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 	Reduceloadx(pf, Kfs, ys); delete Kfs;
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_newton.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_newton.cpp	(revision 16125)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_newton.cpp	(revision 16126)
@@ -56,5 +56,5 @@
 		/*Solver forward model*/
 		if(count==1 || newton==2){
-			femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
+			SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
 			CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 			Reduceloadx(pf,Kfs,ys);delete Kfs;
@@ -68,5 +68,5 @@
 
 		/*Prepare next iteration using Newton's method*/
-		femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);delete df;
+		SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);delete df;
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		Reduceloadx(pf,Kfs,ys);delete Kfs;
@@ -76,5 +76,5 @@
 		pJf->Scale(-1.0); pJf->AXPY(pf,+1.0);
 
-		femmodel->CreateJacobianMatrixx(&Jff,kmax);
+		CreateJacobianMatrixx(&Jff,femmodel,kmax);
 		Solverx(&duf,Jff,pJf,NULL,NULL,femmodel->parameters); delete Jff; delete pJf;
 		uf->AXPY(duf, 1.0); delete duf;
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp	(revision 16125)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp	(revision 16126)
@@ -60,5 +60,5 @@
 		delete ug;
 
-		femmodel->SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL);
+		SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		Reduceloadx(pf, Kfs, ys); delete Kfs;
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp	(revision 16125)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp	(revision 16126)
@@ -61,5 +61,5 @@
 
 		/*solve: */
-		femmodel->SystemMatricesx(&Kff_horiz, &Kfs_horiz, &pf_horiz, &df_horiz, NULL);
+		SystemMatricesx(&Kff_horiz, &Kfs_horiz, &pf_horiz, &df_horiz, NULL,femmodel);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		Reduceloadx(pf_horiz, Kfs_horiz, ys); delete Kfs_horiz;
@@ -75,5 +75,5 @@
 
 		/*solve: */
-		femmodel->SystemMatricesx(&Kff_vert, &Kfs_vert, &pf_vert,  &df_vert,NULL);
+		SystemMatricesx(&Kff_vert, &Kfs_vert, &pf_vert, &df_vert,NULL,femmodel);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		Reduceloadx(pf_vert, Kfs_vert, ys); delete Kfs_vert;
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp	(revision 16125)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp	(revision 16126)
@@ -48,5 +48,5 @@
 
 		delete tf_old; tf_old=tf;
-		femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset);
+		SystemMatricesx(&Kff, &Kfs, &pf,&df, &melting_offset,femmodel);
 		CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
 		Reduceloadx(pf, Kfs, ys); delete Kfs;