Index: ../trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 13880)
+++ ../trunk-jpl/src/c/classes/FemModel.cpp	(revision 13881)
@@ -339,9 +339,12 @@
 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;
+	int  fsize,ssize,flocalsize,slocalsize;
+	int  connectivity, numberofdofspernode;
+	int  analysis_type,configuration_type;
+	int  m,n,M,N;
+	int *d_nnz = NULL;
+	int *o_nnz = NULL;
 
 	/*output*/
 	Matrix<IssmDouble> *Kff  = NULL;
@@ -357,8 +360,11 @@
 	this->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
 
 	/*retrieve node info*/
-	fsize=this->nodes->NumberOfDofs(configuration_type,FsetEnum);
-	ssize=this->nodes->NumberOfDofs(configuration_type,SsetEnum);
+	fsize      = this->nodes->NumberOfDofs(configuration_type,FsetEnum);
+	ssize      = this->nodes->NumberOfDofs(configuration_type,SsetEnum);
+	flocalsize = this->nodes->NumberOfDofsLocal(analysis_type,FsetEnum);
+	slocalsize = this->nodes->NumberOfDofsLocal(analysis_type,SsetEnum);
+
 	numberofdofspernode=this->nodes->MaxNumDofs(configuration_type,GsetEnum);
 
 	if(oldalloc){
@@ -368,7 +374,25 @@
 		pf =new Vector<IssmDouble>(fsize);
 	}
 	else{
-		/*IN PROGRESS*/
+		/*Kff*/
+		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);
+
+		/*Kfs*/
+		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);
+
+		/*Vectors*/
+		df =new Vector<IssmDouble>(flocalsize,fsize);
+		pf =new Vector<IssmDouble>(flocalsize,fsize);
 	}
 
 	/*Allocate output pointers*/
@@ -378,6 +402,136 @@
 	*ppf  = pf;
 }
 /*}}}*/
+void FemModel::MatrixNonzeros(int** pd_nnz,int** po_nnz,int set1enum,int set2enum){/*{{{*/
+
+	/*Intermediary*/
+	int      i,j,k,index,count;
+	int      analysis_type,configuration_type;
+	int      fsize,dim;
+	int      d_nz,o_nz;
+	Element *element = NULL;
+	int     *head    = NULL;
+	int     *next    = NULL;
+
+	/*output*/
+	int *d_nnz = NULL;
+	int *o_nnz = NULL;
+
+	/*retrive parameters: */
+	this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+	this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+	this->parameters->FindParam(&dim,MeshDimensionEnum);
+
+	/*Get vector size and number of nodes*/
+	int numnodes            = nodes->NumberOfNodes(analysis_type);
+	int numberofdofspernode = nodes->MaxNumDofs(configuration_type,GsetEnum);
+	int m                   = nodes->NumberOfDofsLocal(analysis_type,set1enum);
+
+	/*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*/
+	if(dim==2){
+		head=xNew<int>(numnodes); for(i=0;i<numnodes;i++) head[i]=-1;
+		next=xNew<int>(elements->Size()*3); /*3 = number of nodes per element*/
+	}
+	else if(dim==3){
+		head=xNew<int>(numnodes); for(i=0;i<numnodes;i++) head[i]=-1;
+		next=xNew<int>(elements->Size()*6); /*6 = number of nodes per element*/
+	}
+	else{
+		_error_("dim "<<dim<<" not supported yet");
+	}
+	k=0;
+	for(i=0;i<elements->Size();i++){
+		element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+		for(int j=0;j<3;j++){
+			int index =dynamic_cast<Tria*>(element)->nodes[j]->sid;//starts at 0 for a given analysis
+			_assert_(k>=0 && k<numnodes*elements->Size() && index>=0 && index<numnodes);
+			next[k]=head[index];
+			head[index]=k++;
+		}
+	}
+
+	/*OK now count number of dofs and flag each nodes for each node i*/
+	bool *flags                  = xNew<bool>(numnodes);
+	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);
+
+	/*Create connectivity vector*/
+	for(i=0;i<nodes->Size();i++){
+		Node* node=dynamic_cast<Node*>(nodes->GetObjectByOffset(i));
+		if(node->InAnalysis(analysis_type)){
+
+			/*Reinitialize flags to 0*/
+			for(j=0;j<numnodes;j++) flags[j]=false;
+
+			/*Loop over elements that hold node number i*/
+			_assert_(head[node->Sid()]!=-1);
+			for(j=head[node->Sid()];j!=-1;j=next[j]){
+				if(dim==2){
+					index = (int)(double(j)/3);
+				}
+				else if(dim==3){
+					index = (int)(double(j)/6);
+				}
+				else{
+					_error_("dim "<<dim<<" not supported yet");
+				}
+				element=dynamic_cast<Element*>(elements->GetObjectByOffset(index));
+				element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,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>(head);
+	xDelete<int>(next);
+
+	/*sum over all cpus*/
+	MPI_Allreduce((void*)connectivity_clone,(void*)all_connectivity_clone,numnodes,MPI_INT,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(analysis_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]>m)   d_nnz[count]=m;
+					if(o_nnz[count]>fsize-m) o_nnz[count]=fsize-m;
+					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;
+
+}/*}}}*/
 int  FemModel::UpdateVertexPositionsx(void){ /*{{{*/
 
 	int     i;
@@ -627,6 +781,9 @@
 	Kff->Assemble();
 	Kfs->Assemble();
 	pf->Assemble();
+	//Kff->AllocationInfo();
+	//Kfs->AllocationInfo();
+	//_error_("STOP");
 
 	/*Assign output pointers: */
 	if(pKff) *pKff=Kff;
Index: ../trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- ../trunk-jpl/src/c/classes/FemModel.h	(revision 13880)
+++ ../trunk-jpl/src/c/classes/FemModel.h	(revision 13881)
@@ -54,6 +54,7 @@
 		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);
+		void MatrixNonzeros(int** pd_nnz,int** po_nnz,int set1enum,int set2enum);
 		void Solve(void);
 		void OutputResults(void);
 		void SetStaticComm();