Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp	(revision 26130)
+++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp	(revision 26131)
@@ -227,9 +227,9 @@
 	/*Transpose and return Ke*/
 	Ke->Transpose();
-	_assert_(Ke->nrows == numvertices && Ke->ncols == numvertices);
+	_assert_(Ke->nrows == numvertices);
 
 	for(int i=0;i<numvertices;i++){
 		for(int j=0;j<numvertices;j++){
-			ge[i] += Ke->values[i*Ke->ncols + j] * lambda[j];
+			ge[i] += Ke->values[i*Ke->nrows + j] * lambda[j];
 		}
 		//ge[i]=-lambda[i];
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 26130)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 26131)
@@ -947,5 +947,5 @@
 	/*First, figure out size of doflist and create it: */
 	int numberofdofs=0;
-	for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
+	for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,GsetEnum);
 
 	/*Allocate output*/
@@ -955,6 +955,6 @@
 	int count=0;
 	for(int i=0;i<numnodes;i++){
-		nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
-		count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
+		nodes[i]->GetDofList(&doflist[count],approximation_enum,setenum);
+		count+=nodes[i]->GetNumberOfDofs(approximation_enum,GsetEnum);
 	}
 
@@ -970,5 +970,5 @@
 	/*First, figure out size of doflist and create it: */
 	int numberofdofs=0;
-	for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
+	for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
 
 	/*Allocate output*/
@@ -978,6 +978,7 @@
 	int count=0;
 	for(int i=0;i<numnodes;i++){
+		_assert_(count<numberofdofs);
 		nodes[i]->GetDofListLocal(doflist+count,approximation_enum,setenum);
-		count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
+		count+=nodes[i]->GetNumberOfDofs(approximation_enum,GsetEnum);
 	}
 
@@ -4373,6 +4374,6 @@
 
 	/*Copy current stiffness matrix*/
-	values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
-	for(i=0;i<Ke->nrows;i++) for(j=0;j<Ke->ncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j];
+	values=xNew<IssmDouble>(Ke->nrows*Ke->nrows);
+	for(i=0;i<Ke->nrows;i++) for(j=0;j<Ke->nrows;j++) values[i*Ke->nrows+j]=Ke->values[i*Ke->nrows+j];
 
 	/*Get Coordinate Systems transform matrix*/
@@ -4381,5 +4382,5 @@
 	/*Transform matrix: R*Ke*R^T */
 	TripleMultiply(transform,numdofs,numdofs,0,
-				values,Ke->nrows,Ke->ncols,0,
+				values,Ke->nrows,Ke->nrows,0,
 				transform,numdofs,numdofs,1,
 				&Ke->values[0],0);
@@ -4549,6 +4550,6 @@
 
 	/*Copy current stiffness matrix*/
-	values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
-	for(int i=0;i<Ke->nrows*Ke->ncols;i++) values[i]=Ke->values[i];
+	values=xNew<IssmDouble>(Ke->nrows*Ke->nrows);
+	for(int i=0;i<Ke->nrows*Ke->nrows;i++) values[i]=Ke->values[i];
 
 	/*Get Coordinate Systems transform matrix*/
@@ -4557,5 +4558,5 @@
 	/*Transform matrix: R^T*Ke*R */
 	TripleMultiply(transform,numdofs,numdofs,1,
-				values,Ke->nrows,Ke->ncols,0,
+				values,Ke->nrows,Ke->nrows,0,
 				transform,numdofs,numdofs,0,
 				&Ke->values[0],0);
Index: /issm/trunk-jpl/src/c/classes/Node.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Node.cpp	(revision 26130)
+++ /issm/trunk-jpl/src/c/classes/Node.cpp	(revision 26131)
@@ -68,11 +68,11 @@
 		this->svalues        = xNew<IssmDouble>(this->gsize);
 		this->gdoflist       = xNew<int>(this->gsize);
-		this->gdoflist_local = xNewZeroInit<int>(this->gsize);
+		this->gdoflist_local = xNew<int>(this->gsize);
 		this->fsize          = -1;
 		this->ssize          = -1;
-		this->fdoflist       = NULL;
-		this->sdoflist       = NULL;
-		this->fdoflist_local = NULL;
-		this->sdoflist_local = NULL;
+		this->fdoflist       = xNew<int>(this->gsize);
+		this->sdoflist       = xNew<int>(this->gsize);
+		this->fdoflist_local = xNew<int>(this->gsize);
+		this->sdoflist_local = xNew<int>(this->gsize);
 	}
 	else{
@@ -96,4 +96,9 @@
 		this->svalues[i]  = 0.;
 		this->gdoflist[i] = -1;
+		this->fdoflist[i] = -1;
+		this->sdoflist[i] = -1;
+		this->gdoflist_local[i] = -1;
+		this->fdoflist_local[i] = -1;
+		this->sdoflist_local[i] = -1;
 	}
 
@@ -224,10 +229,10 @@
 	}
 	if(output->fsize>0){
-		output->fdoflist=xNew<int>(output->fsize);
-		output->fdoflist_local=xNew<int>(output->fsize);
+		output->fdoflist=xNew<int>(output->gsize);
+		output->fdoflist_local=xNew<int>(output->gsize);
 	}
 	if(output->ssize>0){
-		output->sdoflist=xNew<int>(output->ssize);
-		output->sdoflist_local=xNew<int>(output->ssize);
+		output->sdoflist=xNew<int>(output->gsize);
+		output->sdoflist_local=xNew<int>(output->gsize);
 	}
 
@@ -239,12 +244,8 @@
 		memcpy(output->gdoflist,this->gdoflist,output->gsize*sizeof(int));
 		memcpy(output->gdoflist_local,this->gdoflist_local,output->gsize*sizeof(int));
-	}
-	if(output->fsize>0){
-		memcpy(output->fdoflist,this->fdoflist,output->fsize*sizeof(int));
-		memcpy(output->fdoflist_local,this->fdoflist_local,output->fsize*sizeof(int));
-	}
-	if(output->ssize>0){
-		memcpy(output->sdoflist,this->sdoflist,output->ssize*sizeof(int));
-		memcpy(output->sdoflist_local,this->sdoflist_local,output->ssize*sizeof(int));
+		memcpy(output->fdoflist,this->fdoflist,output->gsize*sizeof(int));
+		memcpy(output->fdoflist_local,this->fdoflist_local,output->gsize*sizeof(int));
+		memcpy(output->sdoflist,this->sdoflist,output->gsize*sizeof(int));
+		memcpy(output->sdoflist_local,this->sdoflist_local,output->gsize*sizeof(int));
 	}
 
@@ -277,9 +278,9 @@
 	marshallhandle->call(this->doftype,gsize);
 	marshallhandle->call(this->gdoflist,gsize);
-	marshallhandle->call(this->fdoflist,fsize);
-	marshallhandle->call(this->sdoflist,ssize);
+	marshallhandle->call(this->fdoflist,gsize);
+	marshallhandle->call(this->sdoflist,gsize);
 	marshallhandle->call(this->gdoflist_local,gsize);
-	marshallhandle->call(this->fdoflist_local,fsize);
-	marshallhandle->call(this->sdoflist_local,ssize);
+	marshallhandle->call(this->fdoflist_local,gsize);
+	marshallhandle->call(this->sdoflist_local,gsize);
 } /*}}}*/
 
@@ -327,16 +328,16 @@
 	_printf_("\n");
 
-	_printf_("   f_doflist (" << this->fsize << "): |");
-	for(i=0;i<this->fsize;i++) _printf_(" " << fdoflist[i] << " |");
+	_printf_("   f_doflist (" << this->gsize << "): |");
+	for(i=0;i<this->gsize;i++) _printf_(" " << fdoflist[i] << " |");
 	_printf_("\n");
-	_printf_("   f_doflist_local (" << this->fsize << "): |");
-	for(i=0;i<this->fsize;i++) _printf_(" " << fdoflist_local[i] << " |");
+	_printf_("   f_doflist_local (" << this->gsize << "): |");
+	for(i=0;i<this->gsize;i++) _printf_(" " << fdoflist_local[i] << " |");
 	_printf_("\n");
 
-	_printf_("   s_doflist (" << this->ssize << "): |");
-	for(i=0;i<this->ssize;i++) _printf_(" " << sdoflist[i] << " |");
+	_printf_("   s_doflist (" << this->gsize << "): |");
+	for(i=0;i<this->gsize;i++) _printf_(" " << sdoflist[i] << " |");
 	_printf_("\n");
-	_printf_("   s_doflist_local (" << this->ssize << "): |");
-	for(i=0;i<this->ssize;i++) _printf_(" " << sdoflist_local[i] << " |");
+	_printf_("   s_doflist_local (" << this->gsize << "): |");
+	for(i=0;i<this->gsize;i++) _printf_(" " << sdoflist_local[i] << " |");
 	_printf_("\n");
 
@@ -386,9 +387,9 @@
 	}
 	else if(setenum==FsetEnum){
-		_assert_(dofindex>=0 && dofindex<fsize);
+		_assert_(dofindex>=0 && dofindex<gsize);
 		return fdoflist[dofindex];
 	}
 	else if(setenum==SsetEnum){
-		_assert_(dofindex>=0 && dofindex<ssize);
+		_assert_(dofindex>=0 && dofindex<gsize);
 		return sdoflist[dofindex];
 	}
@@ -397,226 +398,48 @@
 } /*}}}*/
 void Node::GetDofList(int* outdoflist,int approximation_enum,int setenum){/*{{{*/
+	_assert_(!this->indexingupdate);
 	int i;
-	int count=0;
-	int count2=0;
+
+	int* doflistpointer = NULL;
+	if(setenum==GsetEnum) doflistpointer = gdoflist;
+	else if(setenum==FsetEnum)for(i=0;i<this->gsize;i++) doflistpointer = fdoflist;
+	else if(setenum==SsetEnum)for(i=0;i<this->gsize;i++) doflistpointer = sdoflist;
+	else _error_("not supported");
+
+	if(approximation_enum==NoneApproximationEnum){
+		for(i=0;i<this->gsize;i++) outdoflist[i]=doflistpointer[i];
+	}
+	else{
+		if(doftype){
+			int count = 0;
+			for(i=0;i<this->gsize;i++){
+				if(doftype[i]==approximation_enum) outdoflist[count++]=doflistpointer[i];
+			}
+		}
+		else for(i=0;i<this->gsize;i++) outdoflist[i]=doflistpointer[i];
+	}
+}/*}}}*/
+void Node::GetDofListLocal(int* outdoflist,int approximation_enum,int setenum){/*{{{*/
 
 	_assert_(!this->indexingupdate);
+	int i;
+
+	int* doflistpointer = NULL;
+	if(setenum==GsetEnum) doflistpointer = gdoflist_local;
+	else if(setenum==FsetEnum)for(i=0;i<this->gsize;i++) doflistpointer = fdoflist_local;
+	else if(setenum==SsetEnum)for(i=0;i<this->gsize;i++) doflistpointer = sdoflist_local;
+	else _error_("not supported");
 
 	if(approximation_enum==NoneApproximationEnum){
-		if(setenum==GsetEnum)for(i=0;i<this->gsize;i++) outdoflist[i]=gdoflist[i];
-		if(setenum==FsetEnum)for(i=0;i<this->fsize;i++) outdoflist[i]=fdoflist[i];
-		if(setenum==SsetEnum)for(i=0;i<this->ssize;i++) outdoflist[i]=sdoflist[i];
+		for(i=0;i<this->gsize;i++) outdoflist[i]=doflistpointer[i];
 	}
 	else{
-
-		if(setenum==GsetEnum){
-			if(doftype){
-				count=0;
-				for(i=0;i<this->gsize;i++){
-					if(doftype[i]==approximation_enum){
-						outdoflist[count]=gdoflist[i];
-						count++;
-					}
-				}
-				_assert_(count); //at least one dof should be the approximation requested
-			}
-			else for(i=0;i<this->gsize;i++) outdoflist[i]=gdoflist[i];
-		}
-		else if(setenum==FsetEnum){
-			if(doftype){
-				count=0;
-				count2=0;
-				for(i=0;i<this->gsize;i++){
-					if(f_set[i]){
-						if(doftype[i]==approximation_enum){
-							outdoflist[count]=fdoflist[count2];
-							count++;
-						}
-						count2++;
-					}
-				}
-			}
-			else for(i=0;i<this->fsize;i++) outdoflist[i]=fdoflist[i];
-		}
-		else if(setenum==SsetEnum){
-			if(doftype){
-				count=0;
-				count2=0;
-				for(i=0;i<this->gsize;i++){
-					if(s_set[i]){
-						if(doftype[i]==approximation_enum){
-							outdoflist[count]=sdoflist[count2];
-							count++;
-						}
-						count2++;
-					}
-				}
-			}
-			else for(i=0;i<this->ssize;i++) outdoflist[i]=sdoflist[i];
-		}
-		else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
-	}
-}
-/*}}}*/
-void Node::GetDofListLocal(int* outdoflist,int approximation_enum,int setenum){/*{{{*/
-	int i;
-	int count=0;
-	int count2=0;
-
-	_assert_(!this->indexingupdate);
-
-	if(approximation_enum==NoneApproximationEnum){
-		if(setenum==GsetEnum)for(i=0;i<this->gsize;i++) outdoflist[i]=gdoflist_local[i];
-		if(setenum==FsetEnum)for(i=0;i<this->fsize;i++) outdoflist[i]=fdoflist_local[i];
-		if(setenum==SsetEnum)for(i=0;i<this->ssize;i++) outdoflist[i]=sdoflist_local[i];
-	}
-	else{
-
-		if(setenum==GsetEnum){
-			if(doftype){
-				count=0;
-				for(i=0;i<this->gsize;i++){
-					if(doftype[i]==approximation_enum){
-						outdoflist[count]=gdoflist_local[i];
-						count++;
-					}
-				}
-				_assert_(count); //at least one dof should be the approximation requested
-			}
-			else for(i=0;i<this->gsize;i++) outdoflist[i]=gdoflist_local[i];
-		}
-		else if(setenum==FsetEnum){
-			if(doftype){
-				count=0;
-				count2=0;
-				for(i=0;i<this->gsize;i++){
-					if(f_set[i]){
-						if(doftype[i]==approximation_enum){
-							outdoflist[count]=fdoflist_local[count2];
-							count++;
-						}
-						count2++;
-					}
-				}
-			}
-			else for(i=0;i<this->fsize;i++) outdoflist[i]=fdoflist_local[i];
-		}
-		else if(setenum==SsetEnum){
-			if(doftype){
-				count=0;
-				count2=0;
-				for(i=0;i<this->gsize;i++){
-					if(s_set[i]){
-						if(doftype[i]==approximation_enum){
-							outdoflist[count]=sdoflist_local[count2];
-							count++;
-						}
-						count2++;
-					}
-				}
-			}
-			else for(i=0;i<this->ssize;i++) outdoflist[i]=sdoflist_local[i];
-		}
-		else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
-	}
-}
-/*}}}*/
-void Node::GetLocalDofList(int* outdoflist,int approximation_enum,int setenum){/*{{{*/
-	int i;
-	int count=0;
-	int count2=0;
-
-	_assert_(!this->indexingupdate);
-
-	if(approximation_enum==NoneApproximationEnum){
-		if(setenum==GsetEnum)for(i=0;i<this->gsize;i++) outdoflist[i]=i;
-		else if(setenum==FsetEnum){
-			count=0;
+		if(doftype){
+			int count =0;
 			for(i=0;i<this->gsize;i++){
-				if(f_set[i]){
-					outdoflist[count]=i;
-					count++;
-				}
-			}
-		}
-		else if(setenum==SsetEnum){
-			count=0;
-			for(i=0;i<this->gsize;i++){
-				if(s_set[i]){
-					outdoflist[count]=i;
-					count++;
-				}
-			}
-		}
-		else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
-	}
-	else{
-
-		if(setenum==GsetEnum){
-			if(doftype){
-				count=0;
-				for(i=0;i<this->gsize;i++){
-					if(doftype[i]==approximation_enum){
-						outdoflist[count]=count;
-						count++;
-					}
-				}
-				_assert_(count);
-			}
-			else for(i=0;i<this->gsize;i++) outdoflist[i]=i;
-		}
-		else if(setenum==FsetEnum){
-
-			if(doftype){
-				count=0;
-				count2=0;
-				for(i=0;i<this->gsize;i++){
-					if(doftype[i]==approximation_enum){
-						if(f_set[i]){
-							outdoflist[count]=count2;
-							count++;
-						}
-						count2++;
-					}
-				}
-				_assert_(count2);
-			}
-			else{
-
-				count=0;
-				for(i=0;i<this->gsize;i++){
-					if(f_set[i]){
-						outdoflist[count]=i;
-						count++;
-					}
-				}
-			}
-		}
-		else if(setenum==SsetEnum){
-			if(doftype){
-				count=0;
-				count2=0;
-				for(i=0;i<this->gsize;i++){
-					if(doftype[i]==approximation_enum){
-						if(s_set[i]){
-							outdoflist[count]=count2;
-							count++;
-						}
-						count2++;
-					}
-				}
-				_assert_(count2);
-			}
-			else{
-				count=0;
-				for(i=0;i<this->gsize;i++){
-					if(s_set[i]){
-						outdoflist[count]=i;
-						count++;
-					}
-				}
-			}
-		}
-		else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
+				if(doftype[i]==approximation_enum) outdoflist[count++]=doflistpointer[i];
+			}
+		}
+		else for(i=0;i<this->gsize;i++) outdoflist[i]=doflistpointer[i];
 	}
 }
@@ -660,29 +483,10 @@
 void Node::CreateNodalConstraints(Vector<IssmDouble>* ys){/*{{{*/
 
-	int i;
-	IssmDouble* values=NULL;
-	int count;
-
-	/*Recover values for s set and plug them in constraints vector: */
 	if(this->ssize){
-		values=xNew<IssmDouble>(this->ssize);
-		count=0;
-		for(i=0;i<this->gsize;i++){
-			if(this->s_set[i]){
-				values[count]=this->svalues[i];
-				_assert_(!xIsNan<IssmDouble>(values[count]));
-				count++;
-			}
-		}
-
 		/*Add values into constraint vector: */
-		ys->SetValues(this->ssize,this->sdoflist,values,INS_VAL);
-	}
-
-	/*Free ressources:*/
-	xDelete<IssmDouble>(values);
-
-}
-/*}}}*/
+		ys->SetValues(this->gsize,this->sdoflist,this->svalues,INS_VAL);
+	}
+
+}/*}}}*/
 void Node::Deactivate(void){/*{{{*/
 
@@ -705,9 +509,9 @@
 	_assert_(this->active);
 
-	if(this->f_set[dof] == 0){
+	if(!this->f_set[dof]){
 		if(this->freeze) _error_("Cannot change dof of frozen node");
 		this->indexingupdate = true;
-		this->f_set[dof]=1; 
-		this->s_set[dof]=0;
+		this->f_set[dof]=true; 
+		this->s_set[dof]=false;
 	}
 }
@@ -719,9 +523,9 @@
 	_assert_(dof<this->gsize);
 
-	if(this->f_set[dof] == 1){
+	if(this->f_set[dof]){
 		//if(this->freeze) _error_("Cannot change dof of frozen node");
 		this->indexingupdate = true;
-		this->f_set[dof]=0; //n splits into f (for which we solve) and s (single point constraints)
-		this->s_set[dof]=1;
+		this->f_set[dof]=false; //n splits into f (for which we solve) and s (single point constraints)
+		this->s_set[dof]=true;
 	}
 }
@@ -800,21 +604,12 @@
 /*}}}*/
 bool Node::IsActive(void){/*{{{*/
-
 	return active;
-
-}
-/*}}}*/
+}/*}}}*/
 int  Node::IsClone(){/*{{{*/
-
 	return clone;
-
-}
-/*}}}*/
+}/*}}}*/
 void Node::ReindexingDone(void){/*{{{*/
-
 	this->indexingupdate = false;
-
-}
-/*}}}*/
+}/*}}}*/
 void Node::RelaxConstraint(int dof){/*{{{*/
 
@@ -845,7 +640,7 @@
 			if(this->f_set[i]){
 				_assert_(local_uf);
-				_assert_(this->fdoflist[count]==indices_uf[this->fdoflist_local[count]]);
-
-				values[count] =local_uf[this->fdoflist_local[count]];
+				_assert_(this->fdoflist[i]==indices_uf[this->fdoflist_local[i]]);
+
+				values[count] =local_uf[this->fdoflist_local[i]];
 				indices[count]=this->gdoflist[i];
 				count++;
@@ -865,7 +660,7 @@
 			if(this->s_set[i]){
 				_assert_(local_ys);
-				_assert_(this->sdoflist[count]==indices_ys[this->sdoflist_local[count]]);
-
-				values[count] =local_ys[this->sdoflist_local[count]];
+				_assert_(this->sdoflist[i]==indices_ys[this->sdoflist_local[i]]);
+
+				values[count] =local_ys[this->sdoflist_local[i]];
 				indices[count]=this->gdoflist[i];
 				count++;
@@ -894,10 +689,12 @@
 				_assert_(local_ug);
 				_assert_(this->gdoflist[i]==indices_ug[this->gdoflist_local[i]]);
+				_assert_(this->fdoflist[i]>=0);
 
 				values[count] =local_ug[this->gdoflist_local[i]];
-				indices[count]=this->fdoflist[count];
+				indices[count]=this->fdoflist[i];
 				count++;
 			}
 		}
+		_assert_(count==this->fsize);
 		uf->SetValues(this->fsize,indices,values,INS_VAL);
 
@@ -917,30 +714,8 @@
 		for(int i=0;i<this->gsize;i++) if(f_set[i])size++;
 		this->fsize=size;
-		xDelete<int>(this->fdoflist);
-		xDelete<int>(this->fdoflist_local);
-
-		if(this->fsize){
-			this->fdoflist=xNew<int>(size);
-			this->fdoflist_local=xNew<int>(size);
-		}
-		else{
-			this->fdoflist=NULL;
-			this->fdoflist_local=NULL;
-		}
 	}
 	else if(setenum==SsetEnum){
 		for(int i=0;i<this->gsize;i++) if(s_set[i])size++;
 		this->ssize=size;
-		xDelete<int>(this->sdoflist);
-		xDelete<int>(this->sdoflist_local);
-
-		if(this->ssize){
-			this->sdoflist=xNew<int>(size);
-			this->sdoflist_local=xNew<int>(size);
-		}
-		else{
-			this->sdoflist=NULL;
-			this->sdoflist_local=NULL;
-		}
 	}
 	/*TODO, we should never be here for GSET! add an assert and track down whether this happens*/
@@ -955,16 +730,19 @@
 	if(setenum==GsetEnum){
 		_assert_(this->gsize==0 || this->gdoflist_local);
-		for(int i=0;i<this->gsize;i++) gdoflist_local[i]=dofcount+i;
-		dofcount+=this->gsize;
+		for(int i=0;i<this->gsize;i++) gdoflist_local[i]=dofcount++;
 	}
 	else if(setenum==FsetEnum){
 		_assert_(this->fsize==0 || this->fdoflist_local);
-		for(int i=0;i<this->fsize;i++) fdoflist_local[i]=dofcount+i;
-		dofcount+=this->fsize;
+		for(int i=0;i<this->gsize;i++){
+			if(this->f_set[i]) fdoflist_local[i]=dofcount++;
+			else               fdoflist_local[i]=-1;
+		}
 	}
 	else if(setenum==SsetEnum){
 		_assert_(this->ssize==0 || this->sdoflist_local);
-		for(int i=0;i<this->ssize;i++) sdoflist_local[i]=dofcount+i;
-		dofcount+=this->ssize;
+		for(int i=0;i<this->gsize;i++){
+			if(this->s_set[i]) sdoflist_local[i]=dofcount++;
+			else               sdoflist_local[i]=-1;
+		}
 	}
 	else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
@@ -972,6 +750,5 @@
 	/*Assign output pointers: */
 	*pdofcount=dofcount;
-}
-/*}}}*/
+}/*}}}*/
 void Node::DistributeGlobalDofsMasters(int dofcount,int setenum){/*{{{*/
 
@@ -986,9 +763,15 @@
 	else if(setenum==FsetEnum){
 		_assert_(this->fsize==0 || this->fdoflist);
-		for(int i=0;i<this->fsize;i++) this->fdoflist[i]=this->fdoflist_local[i]+dofcount;
+		for(int i=0;i<this->gsize;i++){
+			if(this->f_set[i]) fdoflist[i]=this->fdoflist_local[i]+dofcount;
+			else               fdoflist[i]=-1;
+		}
 	}
 	else if(setenum==SsetEnum){
 		_assert_(this->ssize==0 || this->sdoflist);
-		for(int i=0;i<this->ssize;i++) this->sdoflist[i]=this->sdoflist_local[i]+dofcount;
+		for(int i=0;i<this->gsize;i++){
+			if(this->s_set[i]) sdoflist[i]=this->sdoflist_local[i]+dofcount;
+			else               sdoflist[i]=-1;
+		}
 	}
 	else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
@@ -1005,8 +788,8 @@
 			break;
 		case FsetEnum:
-			for(int j=0;j<this->fsize;j++) truedofs[j]=fdoflist[j];
+			for(int j=0;j<this->gsize;j++) truedofs[j]=fdoflist[j];
 			break;
 		case SsetEnum:
-			for(int j=0;j<this->ssize;j++) truedofs[j]=sdoflist[j];
+			for(int j=0;j<this->gsize;j++) truedofs[j]=sdoflist[j];
 			break;
 		default:
@@ -1027,8 +810,8 @@
 			break;
 		case FsetEnum:
-			for(int j=0;j<this->fsize;j++) fdoflist[j]=alltruedofs[j];
+			for(int j=0;j<this->gsize;j++) fdoflist[j]=alltruedofs[j];
 			break;
 		case SsetEnum:
-			for(int j=0;j<this->ssize;j++) sdoflist[j]=alltruedofs[j];
+			for(int j=0;j<this->gsize;j++) sdoflist[j]=alltruedofs[j];
 			break;
 		default:
@@ -1041,6 +824,5 @@
 int* GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation){/*{{{*/
 
-	int  i,numdof,count;
-	int* ndof_list=NULL;
+	/*output*/
 	int *doflist = NULL;
 
@@ -1048,10 +830,10 @@
 
 		/*Allocate:*/
-		ndof_list=xNew<int>(numnodes);
+		int* ndof_list=xNew<int>(numnodes);
 
 		/*First, figure out size of doflist: */
-		numdof=0;
-		for(i=0;i<numnodes;i++){
-			ndof_list[i]=nodes[i]->GetNumberOfDofs(approximation,setenum);
+		int numdof=0;
+		for(int i=0;i<numnodes;i++){
+			ndof_list[i]=nodes[i]->GetNumberOfDofs(approximation,GsetEnum);
 			numdof+=ndof_list[i];
 		}
@@ -1062,6 +844,6 @@
 
 			/*Populate: */
-			count=0;
-			for(i=0;i<numnodes;i++){
+			int count=0;
+			for(int i=0;i<numnodes;i++){
 				nodes[i]->GetDofList(&doflist[count],approximation,setenum);
 				count+=ndof_list[i];
@@ -1069,65 +851,9 @@
 		}
 		else doflist=NULL;
-	}
-	/*Free ressources:*/
-	xDelete<int>(ndof_list);
-
-	return doflist;
-}
-/*}}}*/
-int* GetLocalDofList(Node** nodes,int numnodes,int setenum,int approximation){ /*{{{*/
-
-	int  i,j,count,numdof,numgdof;
-	int* ndof_list=NULL;
-	int* ngdof_list_cumulative=NULL;
-	int *doflist = NULL;
-
-	if(numnodes){
-		/*allocate: */
-		ndof_list=xNew<int>(numnodes);
-		ngdof_list_cumulative=xNew<int>(numnodes);
-
-		/*Get number of dofs per node, and total for this given set*/
-		numdof=0;
-		numgdof=0;
-		for(i=0;i<numnodes;i++){
-
-			/*Cumulative list= number of dofs before node i*/
-			ngdof_list_cumulative[i]=numgdof;
-
-			/*Number of dofs for node i for given set and for the g set*/
-			ndof_list[i]=nodes[i]->GetNumberOfDofs(approximation,setenum);
-			numgdof    +=nodes[i]->GetNumberOfDofs(approximation,GsetEnum);
-			numdof     +=ndof_list[i];
-		}
-
-		if(numdof){
-			/*Allocate: */
-			doflist=xNew<int>(numdof);
-
-			/*Populate: */
-			count=0;
-			for(i=0;i<numnodes;i++){
-				nodes[i]->GetLocalDofList(&doflist[count],approximation,setenum);
-				count+=ndof_list[i];
-			}
-
-			/*We now have something like: [0 1 0 2 1 2]. Offset by gsize, to get something like: [0 1 2 4 6 7]:*/
-			count=0;
-			for(i=0;i<numnodes;i++){
-				for(j=0;j<ndof_list[i];j++){
-					doflist[count+j]+=ngdof_list_cumulative[i];
-				}
-				count+=ndof_list[i];
-			}
-		}
-		else doflist=NULL;
-	}
-
-	/*Free ressources:*/
-	xDelete<int>(ndof_list);
-	xDelete<int>(ngdof_list_cumulative);
-
-	/*CLean-up and return*/
+
+		/*Free ressources:*/
+		xDelete<int>(ndof_list);
+	}
+
 	return doflist;
 }
Index: /issm/trunk-jpl/src/c/classes/Node.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Node.h	(revision 26130)
+++ /issm/trunk-jpl/src/c/classes/Node.h	(revision 26131)
@@ -94,5 +94,4 @@
 		void  GetDofList(int* poutdoflist,int approximation_enum,int setenum);
 		void  GetDofListLocal(int* poutdoflist,int approximation_enum,int setenum);
-		void  GetLocalDofList(int* poutdoflist,int approximation_enum,int setenum);
 		int   GetNumberOfDofs(int approximation_enum,int setenum);
 		void  HardDeactivate(void);
@@ -118,5 +117,4 @@
 /*Methods inherent to Node: */
 int* GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation);
-int* GetLocalDofList(Node** nodes,int numnodes,int setenum,int approximation);
 int  GetNumberOfDofs(Node** nodes,int numnodes,int setenum,int approximation);
 
Index: /issm/trunk-jpl/src/c/classes/Nodes.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Nodes.cpp	(revision 26130)
+++ /issm/trunk-jpl/src/c/classes/Nodes.cpp	(revision 26131)
@@ -187,5 +187,5 @@
 	 * object that is not a clone, tell them to show their dofs, so that later on, they can get picked
 	 * up by their clones: */
-	int  maxdofspernode = this->MaxNumDofs(setenum);
+	int  maxdofspernode = this->MaxNumDofs(GsetEnum);
 	int* truedofs       = xNewZeroInit<int>(this->Size()*maxdofspernode); //only one alloc
 	for(int rank=0;rank<num_procs;rank++){
Index: /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp	(revision 26130)
+++ /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp	(revision 26131)
@@ -21,22 +21,10 @@
 
 	this->nrows=0;
-	this->ncols=0;
 	this->values=NULL;
-	this->dofsymmetrical=false;
-
-	this->row_fsize=0;
-	this->row_flocaldoflist=NULL;
-	this->row_fglobaldoflist=NULL;
-	this->row_ssize=0;
-	this->row_slocaldoflist=NULL;
-	this->row_sglobaldoflist=NULL;
-
-	this->col_fsize=0;
-	this->col_flocaldoflist=NULL;
-	this->col_fglobaldoflist=NULL;
-	this->col_ssize=0;
-	this->col_slocaldoflist=NULL;
-	this->col_sglobaldoflist=NULL;
-
+
+	this->fsize=0;
+	this->fglobaldoflist=NULL;
+	this->ssize=0;
+	this->sglobaldoflist=NULL;
 }
 /*}}}*/
@@ -53,6 +41,4 @@
 	int i,j,counter;
 	int gsize,fsize,ssize;
-	int* P=NULL;
-	bool found;
 
 	/*If one of the two matrix is NULL, we copy the other one*/
@@ -70,13 +56,12 @@
 
 	/*General Case: Ke1 and Ke2 are not empty*/
-	if(!Ke1->dofsymmetrical || !Ke2->dofsymmetrical) _error_("merging 2 non dofsymmetrical matrices not implemented yet");
 
 	/*Initialize itransformation matrix Ke[P[i]] = Ke2[i]*/
-	P=xNew<int>(Ke2->nrows);
+	int* P=xNew<int>(Ke2->nrows);
 
 	/*1: Get the new numbering of Ke2 and get size of the new matrix*/
 	gsize=Ke1->nrows;
 	for(i=0;i<Ke2->nrows;i++){
-		found=false;
+		bool found=false;
 		for(j=0;j<Ke1->nrows;j++){
 			if(Ke2->gglobaldoflist[i]==Ke1->gglobaldoflist[j]){
@@ -91,84 +76,34 @@
 	/*2: Initialize static fields*/
 	this->nrows=gsize;
-	this->ncols=gsize;
-	this->dofsymmetrical=true;
 
 	/*Gset and values*/
 	this->gglobaldoflist=xNew<int>(this->nrows);
-	this->values=xNewZeroInit<IssmDouble>(this->nrows*this->ncols);
+	this->fglobaldoflist=xNew<int>(this->nrows);
+	this->sglobaldoflist=xNew<int>(this->nrows);
+	this->values=xNewZeroInit<IssmDouble>(this->nrows*this->nrows);
 	for(i=0;i<Ke1->nrows;i++){
-		for(j=0;j<Ke1->ncols;j++){
-			this->values[i*this->ncols+j] += Ke1->values[i*Ke1->ncols+j];
+		for(j=0;j<Ke1->nrows;j++){
+			this->values[i*this->nrows+j] += Ke1->values[i*Ke1->nrows+j];
 		}
 		this->gglobaldoflist[i]=Ke1->gglobaldoflist[i];
+		this->fglobaldoflist[i]=Ke1->fglobaldoflist[i];
+		this->sglobaldoflist[i]=Ke1->sglobaldoflist[i];
 	}
 	for(i=0;i<Ke2->nrows;i++){
-		for(j=0;j<Ke2->ncols;j++){
-			this->values[P[i]*this->ncols+P[j]] += Ke2->values[i*Ke2->ncols+j];
+		for(j=0;j<Ke2->nrows;j++){
+			this->values[P[i]*this->nrows+P[j]] += Ke2->values[i*Ke2->nrows+j];
 		}
 		this->gglobaldoflist[P[i]]=Ke2->gglobaldoflist[i];
+		this->fglobaldoflist[P[i]]=Ke2->fglobaldoflist[i];
+		this->sglobaldoflist[P[i]]=Ke2->sglobaldoflist[i];
 	}
 
 	/*Fset*/
-	fsize=Ke1->row_fsize;
-	for(i=0;i<Ke2->row_fsize;i++){
-		if(P[Ke2->row_flocaldoflist[i]] >= Ke1->nrows) fsize++;
-	}
-	this->row_fsize=fsize;
-	if(fsize){
-		this->row_flocaldoflist =xNew<int>(fsize);
-		this->row_fglobaldoflist=xNew<int>(fsize);
-		for(i=0;i<Ke1->row_fsize;i++){
-			this->row_flocaldoflist[i] =Ke1->row_flocaldoflist[i];
-			this->row_fglobaldoflist[i]=Ke1->row_fglobaldoflist[i];
-		}
-		counter=Ke1->row_fsize;
-		for(i=0;i<Ke2->row_fsize;i++){
-			if(P[Ke2->row_flocaldoflist[i]] >= Ke1->nrows){
-				this->row_flocaldoflist[counter] =P[Ke2->row_flocaldoflist[i]];
-				this->row_fglobaldoflist[counter]=Ke2->row_fglobaldoflist[i];
-				counter++;
-			}
-		}
-	}
-	else{
-		this->row_flocaldoflist=NULL;
-		this->row_fglobaldoflist=NULL;
-	}
+	this->fsize=0;
+	for(i=0;i<this->nrows;i++) if(this->fglobaldoflist[i]>=0) this->fsize++;
 
 	/*Sset*/
-	ssize=Ke1->row_ssize;
-	for(i=0;i<Ke2->row_ssize;i++){
-		if(P[Ke2->row_slocaldoflist[i]] >= Ke1->nrows) ssize++;
-	}
-	this->row_ssize=ssize;
-	if(ssize){
-		this->row_slocaldoflist =xNew<int>(ssize);
-		this->row_sglobaldoflist=xNew<int>(ssize);
-		for(i=0;i<Ke1->row_ssize;i++){
-			this->row_slocaldoflist[i] =Ke1->row_slocaldoflist[i];
-			this->row_sglobaldoflist[i]=Ke1->row_sglobaldoflist[i];
-		}
-		counter=Ke1->row_ssize;
-		for(i=0;i<Ke2->row_ssize;i++){
-			if(P[Ke2->row_slocaldoflist[i]] >= Ke1->nrows){
-				this->row_slocaldoflist[counter] =P[Ke2->row_slocaldoflist[i]];
-				this->row_sglobaldoflist[counter]=Ke2->row_sglobaldoflist[i];
-				counter++;
-			}
-		}
-	}
-	else{
-		this->row_slocaldoflist=NULL;
-		this->row_sglobaldoflist=NULL;
-	}
-
-	/*don't do cols, we can pick them up from the rows: */
-	this->col_fsize=0;
-	this->col_flocaldoflist=NULL;
-	this->col_fglobaldoflist=NULL;
-	this->col_ssize=0;
-	this->col_slocaldoflist=NULL;
-	this->col_sglobaldoflist=NULL;
+	this->ssize=0;
+	for(i=0;i<this->nrows;i++) if(this->sglobaldoflist[i]>=0) this->ssize++;
 
 	/*clean-up*/
@@ -193,10 +128,8 @@
 
 	/*get Matrix size and properties*/
-	this->dofsymmetrical=true;
 	this->nrows=GetNumberOfDofs(nodes,numnodes,GsetEnum,approximation);
-	this->ncols=this->nrows;
 
 	/*fill values with 0: */
-	this->values=xNewZeroInit<IssmDouble>(this->nrows*this->ncols);
+	this->values=xNewZeroInit<IssmDouble>(this->nrows*this->nrows);
 
 	/*g list*/
@@ -204,18 +137,8 @@
 
 	/*get dof lists for f and s set: */
-	this->row_fsize=GetNumberOfDofs(nodes,numnodes,FsetEnum,approximation);
-	this->row_flocaldoflist =GetLocalDofList( nodes,numnodes,FsetEnum,approximation);
-	this->row_fglobaldoflist=GetGlobalDofList(nodes,numnodes,FsetEnum,approximation);
-	this->row_ssize=GetNumberOfDofs(nodes,numnodes,SsetEnum,approximation);
-	this->row_slocaldoflist =GetLocalDofList( nodes,numnodes,SsetEnum,approximation);
-	this->row_sglobaldoflist=GetGlobalDofList(nodes,numnodes,SsetEnum,approximation);
-
-	/*Because this matrix is "dofsymmetrical" don't do cols, we can pick them up from the rows: */
-	this->col_fsize=0;
-	this->col_flocaldoflist=NULL;
-	this->col_fglobaldoflist=NULL;
-	this->col_ssize=0;
-	this->col_slocaldoflist=NULL;
-	this->col_sglobaldoflist=NULL;
+	this->fsize=GetNumberOfDofs(nodes,numnodes,FsetEnum,approximation);
+	this->fglobaldoflist=GetGlobalDofList(nodes,numnodes,FsetEnum,approximation);
+	this->ssize=GetNumberOfDofs(nodes,numnodes,SsetEnum,approximation);
+	this->sglobaldoflist=GetGlobalDofList(nodes,numnodes,SsetEnum,approximation);
 }
 /*}}}*/
@@ -224,12 +147,6 @@
 	xDelete<IssmDouble>(this->values);
 	xDelete<int>(this->gglobaldoflist);
-	xDelete<int>(this->row_flocaldoflist);
-	xDelete<int>(this->row_fglobaldoflist);
-	xDelete<int>(this->row_slocaldoflist);
-	xDelete<int>(this->row_sglobaldoflist);
-	xDelete<int>(this->col_flocaldoflist);
-	xDelete<int>(this->col_fglobaldoflist);
-	xDelete<int>(this->col_slocaldoflist);
-	xDelete<int>(this->col_sglobaldoflist);
+	xDelete<int>(this->fglobaldoflist);
+	xDelete<int>(this->sglobaldoflist);
 }
 /*}}}*/
@@ -246,23 +163,16 @@
 	this->CheckConsistency();
 
-	if(this->dofsymmetrical){
-		/*only use row dofs to add values into global matrices: */
-
-		if(this->row_fsize){
-			/*first, retrieve values that are in the f-set from the g-set values matrix: */
-			localvalues=xNew<IssmDouble>(this->row_fsize);
-			for(int i=0;i<this->row_fsize;i++){
-				localvalues[i] = this->values[this->ncols*this->row_flocaldoflist[i]+ this->row_flocaldoflist[i]];
-			}
-
-			/*add local values into global  matrix, using the fglobaldoflist: */
-			pf->SetValues(this->row_fsize,this->row_fglobaldoflist,localvalues,ADD_VAL);
-
-			/*Free ressources:*/
-			xDelete<IssmDouble>(localvalues);
-		}
-	}
-	else{
-		_error_("non dofsymmetrical matrix AddToGlobal routine not support yet!");
+	if(this->fsize){
+		/*first, retrieve values that are in the f-set from the g-set values matrix: */
+		localvalues=xNew<IssmDouble>(this->nrows);
+		for(int i=0;i<this->nrows;i++){
+			localvalues[i] = this->values[this->nrows*i + i];
+		}
+
+		/*add local values into global  matrix, using the fglobaldoflist: */
+		pf->SetValues(this->fsize,this->fglobaldoflist,localvalues,ADD_VAL);
+
+		/*Free ressources:*/
+		xDelete<IssmDouble>(localvalues);
 	}
 
@@ -283,42 +193,13 @@
 	this->CheckConsistency();
 
-	if(this->dofsymmetrical){
-		/*only use row dofs to add values into global matrices: */
-
-		if(this->row_fsize){
-			/*first, retrieve values that are in the f-set from the g-set values matrix: */
-			IssmDouble* localvalues=xNew<IssmDouble>(this->row_fsize*this->row_fsize);
-			for(int i=0;i<this->row_fsize;i++){
-				for(int j=0;j<this->row_fsize;j++){
-					localvalues[this->row_fsize*i+j]=this->values[this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]];
-				}
-			}
-
-			/*add local values into global  matrix, using the fglobaldoflist: */
-			Kff->SetValues(this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,localvalues,ADD_VAL);
-
-			/*Free ressources:*/
-			xDelete<IssmDouble>(localvalues);
-		}
-
-		if((this->row_ssize!=0) && (this->row_fsize!=0)){
-			/*first, retrieve values that are in the f and s-set from the g-set values matrix: */
-			IssmDouble* localvalues=xNew<IssmDouble>(this->row_fsize*this->row_ssize);
-			for(int i=0;i<this->row_fsize;i++){
-				for(int j=0;j<this->row_ssize;j++){
-					localvalues[this->row_ssize*i+j]=this->values[this->ncols*this->row_flocaldoflist[i]+this->row_slocaldoflist[j]];
-				}
-			}
-			/*add local values into global  matrix, using the fglobaldoflist: */
-			Kfs->SetValues(this->row_fsize,this->row_fglobaldoflist,this->row_ssize,this->row_sglobaldoflist,localvalues,ADD_VAL);
-
-			/*Free ressources:*/
-			xDelete<IssmDouble>(localvalues);
-		}
-	}
-	else{
-		_error_("non dofsymmetrical matrix AddToGlobal routine not support yet!");
-	}
-
+	/*only use row dofs to add values into global matrices: */
+	if(this->fsize){
+		/*add local values into global  matrix, using the fglobaldoflist: */
+		Kff->SetValues(this->nrows,this->fglobaldoflist,this->nrows,this->fglobaldoflist,this->values,ADD_VAL);
+	}
+
+	if((this->ssize!=0) && (this->fsize!=0)){
+		Kfs->SetValues(this->nrows,this->fglobaldoflist,this->nrows,this->sglobaldoflist,this->values,ADD_VAL);
+	}
 }
 /*}}}*/
@@ -331,27 +212,7 @@
 	this->CheckConsistency();
 
-	if(this->dofsymmetrical){
-		/*only use row dofs to add values into global matrices: */
-
-		if(this->row_fsize){
-			/*first, retrieve values that are in the f-set from the g-set values matrix: */
-			IssmDouble* localvalues=xNew<IssmDouble>(this->row_fsize*this->row_fsize);
-			for(int i=0;i<this->row_fsize;i++){
-				for(int j=0;j<this->row_fsize;j++){
-					localvalues[this->row_fsize*i+j]=this->values[this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]];
-				}
-			}
-			/*add local values into global  matrix, using the fglobaldoflist: */
-			Jff->SetValues(this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,localvalues,ADD_VAL);
-
-			/*Free ressources:*/
-			xDelete<IssmDouble>(localvalues);
-		}
-
-	}
-	else{
-		_error_("non dofsymmetrical matrix AddToGlobal routine not support yet!");
-	}
-
+	if(this->fsize){
+		Jff->SetValues(this->nrows,this->fglobaldoflist,this->nrows,this->fglobaldoflist,this->values,ADD_VAL);
+	}
 }
 /*}}}*/
@@ -360,8 +221,8 @@
 	#ifdef _ISSM_DEBUG_ 
 	for (int i=0;i<this->nrows;i++){
-		for(int j=0;j<this->ncols;j++){
-			if (xIsNan<IssmDouble>(this->values[i*this->ncols+j])) _error_("NaN found in Element Matrix");
-			if (xIsInf<IssmDouble>(this->values[i*this->ncols+j])) _error_("Inf found in Element Matrix");
-			if (fabs(this->values[i*this->ncols+j])>1.e+50) _error_("Element Matrix values exceeds 1.e+50");
+		for(int j=0;j<this->nrows;j++){
+			if (xIsNan<IssmDouble>(this->values[i*this->nrows+j])) _error_("NaN found in Element Matrix");
+			if (xIsInf<IssmDouble>(this->values[i*this->nrows+j])) _error_("Inf found in Element Matrix");
+			if (fabs(this->values[i*this->nrows+j])>1.e+50) _error_("Element Matrix values exceeds 1.e+50");
 		}
 	}
@@ -374,11 +235,9 @@
 	_printf_("Element Matrix echo:\n");
 	_printf_("   nrows: " << this->nrows << "\n");
-	_printf_("   ncols: " << this->ncols << "\n");
-	_printf_("   dofsymmetrical: " << (dofsymmetrical?"true":"false") << "\n");
 
 	_printf_("   values:\n");
 	for(i=0;i<nrows;i++){
 		_printf_(setw(4) << right << i << ": ");
-		for(j=0;j<ncols;j++) _printf_( " " << setw(11) << setprecision (5) << right << values[i*ncols+j]);
+		for(j=0;j<nrows;j++) _printf_( " " << setw(11) << setprecision (5) << right << values[i*nrows+j]);
 		_printf_("\n");
 	}
@@ -387,29 +246,11 @@
 	if(gglobaldoflist) for(i=0;i<nrows;i++) _printf_(" " << gglobaldoflist[i]); _printf_("\n");
 
-	_printf_("   row_fsize: " << row_fsize << "\n");
-	_printf_("   row_flocaldoflist  (" << row_flocaldoflist << "): ");
-	if(row_flocaldoflist) for(i=0;i<row_fsize;i++) _printf_(" " << row_flocaldoflist[i] ); _printf_(" \n");
-	_printf_("   row_fglobaldoflist  (" << row_fglobaldoflist << "): ");
-	if(row_fglobaldoflist) for(i=0;i<row_fsize;i++)_printf_(" " << row_fglobaldoflist[i]); _printf_(" \n");
-
-	_printf_("   row_ssize: " << row_ssize << "\n");
-	_printf_("   row_slocaldoflist  (" << row_slocaldoflist << "): ");
-	if(row_slocaldoflist) for(i=0;i<row_ssize;i++) _printf_(" " << row_slocaldoflist[i] ); _printf_(" \n");
-	_printf_("   row_sglobaldoflist  (" << row_sglobaldoflist << "): ");
-	if(row_sglobaldoflist) for(i=0;i<row_ssize;i++)_printf_(" " << row_sglobaldoflist[i]); _printf_(" \n");
-
-	if(!dofsymmetrical){
-		_printf_("   col_fsize: " << col_fsize << "\n");
-		_printf_("   col_flocaldoflist  (" << col_flocaldoflist << "): ");
-		if(col_flocaldoflist) for(i=0;i<col_fsize;i++) _printf_(" " << col_flocaldoflist[i] ); _printf_(" \n");
-		_printf_("   col_fglobaldoflist  (" << col_fglobaldoflist << "): ");
-		if(col_fglobaldoflist) for(i=0;i<col_fsize;i++)_printf_(" " << col_fglobaldoflist[i]); _printf_(" \n");
-
-		_printf_("   col_ssize: " << col_ssize << "\n");
-		_printf_("   col_slocaldoflist  (" << col_slocaldoflist << "): ");
-		if(col_slocaldoflist) for(i=0;i<col_ssize;i++) _printf_(" " << col_slocaldoflist[i] ); _printf_(" \n");
-		_printf_("   col_sglobaldoflist  (" << col_sglobaldoflist << "): ");
-		if(col_sglobaldoflist) for(i=0;i<col_ssize;i++)_printf_(" " << col_sglobaldoflist[i]); _printf_(" \n");
-	}
+	_printf_("   fsize: " << fsize << "\n");
+	_printf_("   fglobaldoflist  (" << fglobaldoflist << "): ");
+	for(i=0;i<nrows;i++)_printf_(" " << fglobaldoflist[i]); _printf_(" \n");
+
+	_printf_("   ssize: " << ssize << "\n");
+	_printf_("   sglobaldoflist  (" << sglobaldoflist << "): ");
+	for(i=0;i<nrows;i++)_printf_(" " << sglobaldoflist[i]); _printf_(" \n");
 }
 /*}}}*/
@@ -420,71 +261,27 @@
 
 	this->nrows =Ke->nrows;
-	this->ncols =Ke->ncols;
-	this->dofsymmetrical=Ke->dofsymmetrical;
-
-	this->values=xNew<IssmDouble>(this->nrows*this->ncols);
-	xMemCpy<IssmDouble>(this->values,Ke->values,this->nrows*this->ncols);
+
+	this->values=xNew<IssmDouble>(this->nrows*this->nrows);
+	xMemCpy<IssmDouble>(this->values,Ke->values,this->nrows*this->nrows);
 
 	this->gglobaldoflist=xNew<int>(this->nrows);
 	xMemCpy<int>(this->gglobaldoflist,Ke->gglobaldoflist,this->nrows);
 
-	this->row_fsize=Ke->row_fsize;
-	if(this->row_fsize){
-		this->row_flocaldoflist=xNew<int>(this->row_fsize);
-		xMemCpy<int>(this->row_flocaldoflist,Ke->row_flocaldoflist,this->row_fsize);
-		this->row_fglobaldoflist=xNew<int>(this->row_fsize);
-		xMemCpy<int>(this->row_fglobaldoflist,Ke->row_fglobaldoflist,this->row_fsize);
-	}
-	else{
-		this->row_flocaldoflist=NULL;
-		this->row_fglobaldoflist=NULL;
-	}
-
-	this->row_ssize=Ke->row_ssize;
-	if(this->row_ssize){
-		this->row_slocaldoflist=xNew<int>(this->row_ssize);
-		xMemCpy<int>(this->row_slocaldoflist,Ke->row_slocaldoflist,this->row_ssize);
-		this->row_sglobaldoflist=xNew<int>(this->row_ssize);
-		xMemCpy<int>(this->row_sglobaldoflist,Ke->row_sglobaldoflist,this->row_ssize);
-	}
-	else{
-		this->row_slocaldoflist=NULL;
-		this->row_sglobaldoflist=NULL;
-	}
-
-	this->col_fsize=Ke->col_fsize;
-	if(this->col_fsize){
-		this->col_flocaldoflist=xNew<int>(this->col_fsize);
-		xMemCpy<int>(this->col_flocaldoflist,Ke->col_flocaldoflist,this->col_fsize);
-		this->col_fglobaldoflist=xNew<int>(this->col_fsize);
-		xMemCpy<int>(this->col_fglobaldoflist,Ke->col_fglobaldoflist,this->col_fsize);
-	}
-	else{
-		this->col_flocaldoflist=NULL;
-		this->col_fglobaldoflist=NULL;
-	}
-
-	this->col_ssize=Ke->col_ssize;
-	if(this->col_ssize){
-		this->col_slocaldoflist=xNew<int>(this->col_ssize);
-		xMemCpy<int>(this->col_slocaldoflist,Ke->col_slocaldoflist,this->col_ssize);
-		this->col_sglobaldoflist=xNew<int>(this->col_ssize);
-		xMemCpy<int>(this->col_sglobaldoflist,Ke->col_sglobaldoflist,this->col_ssize);
-	}
-	else{
-		this->col_slocaldoflist=NULL;
-		this->col_sglobaldoflist=NULL;
-	}
+	this->fsize=Ke->fsize;
+	this->fglobaldoflist=xNew<int>(this->nrows);
+	xMemCpy<int>(this->fglobaldoflist,Ke->fglobaldoflist,this->nrows);
+
+	this->ssize=Ke->ssize;
+	this->sglobaldoflist=xNew<int>(this->nrows);
+	xMemCpy<int>(this->sglobaldoflist,Ke->sglobaldoflist,this->nrows);
 }
 /*}}}*/
 void ElementMatrix::Lump(void){/*{{{*/
 
-	if(!dofsymmetrical) _error_("not supported yet");
-
 	for(int i=0;i<this->nrows;i++){
-		for(int j=0;j<this->ncols;j++){
+		for(int j=0;j<this->nrows;j++){
 			if(i!=j){
-				this->values[i*this->ncols+i] += this->values[i*this->ncols+j];
-				this->values[i*this->ncols+j]  = 0.;
+				this->values[i*this->nrows+i] += this->values[i*this->nrows+j];
+				this->values[i*this->nrows+j]  = 0.;
 			}
 		}
@@ -499,15 +296,6 @@
 	ElementMatrix* Ke_copy=new ElementMatrix(this);
 
-	/*Update sizes*/
-	this->nrows=Ke_copy->ncols;
-	this->ncols=Ke_copy->nrows;
-
 	/*Transpose values*/
-	for (int i=0;i<this->nrows;i++) for(int j=0;j<this->ncols;j++) this->values[i*this->ncols+j]=Ke_copy->values[j*Ke_copy->ncols+i];
-
-	/*Transpose indices*/
-	if(!dofsymmetrical){
-		_error_("not supported yet");
-	}
+	for (int i=0;i<this->nrows;i++) for(int j=0;j<this->nrows;j++) this->values[i*this->nrows+j]=Ke_copy->values[j*Ke_copy->nrows+i];
 
 	/*Clean up and return*/
@@ -534,5 +322,5 @@
 
 	/*Checks in debugging mode*/
-	_assert_(this->nrows==this->ncols && bsize>0 && bsize<this->ncols && this->values); 
+	_assert_(bsize>0 && bsize<this->nrows && this->values); 
 
 	/*Intermediaries*/
@@ -571,8 +359,8 @@
 	Kbi = xNew<IssmDouble>(bsize*isize);
 	Kbb = xNew<IssmDouble>(bsize*bsize);
-	for(i=0;i<isize;i++) for(j=0;j<isize;j++) Kii[i*isize+j] = this->values[iindices[i]*this->ncols + iindices[j]];
-	for(i=0;i<isize;i++) for(j=0;j<bsize;j++) Kib[i*bsize+j] = this->values[iindices[i]*this->ncols + bindices[j]];
-	for(i=0;i<bsize;i++) for(j=0;j<isize;j++) Kbi[i*isize+j] = this->values[bindices[i]*this->ncols + iindices[j]];
-	for(i=0;i<bsize;i++) for(j=0;j<bsize;j++) Kbb[i*bsize+j] = this->values[bindices[i]*this->ncols + bindices[j]];
+	for(i=0;i<isize;i++) for(j=0;j<isize;j++) Kii[i*isize+j] = this->values[iindices[i]*this->nrows + iindices[j]];
+	for(i=0;i<isize;i++) for(j=0;j<bsize;j++) Kib[i*bsize+j] = this->values[iindices[i]*this->nrows + bindices[j]];
+	for(i=0;i<bsize;i++) for(j=0;j<isize;j++) Kbi[i*isize+j] = this->values[bindices[i]*this->nrows + iindices[j]];
+	for(i=0;i<bsize;i++) for(j=0;j<bsize;j++) Kbb[i*bsize+j] = this->values[bindices[i]*this->nrows + bindices[j]];
 
 	/*Invert Kbb*/
@@ -601,8 +389,8 @@
 
 	/*Update matrix values*/
-	for(i=0;i<this->nrows*this->ncols;i++) this->values[i]=0.;
+	for(i=0;i<this->nrows*this->nrows;i++) this->values[i]=0.;
 	for(i=0;i<isize;i++){
 		for(j=0;j<isize;j++){
-			this->values[iindices[i]*this->ncols + iindices[j]] = Ktemp[i*isize+j];
+			this->values[iindices[i]*this->nrows + iindices[j]] = Ktemp[i*isize+j];
 		}
 	}
Index: /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.h	(revision 26130)
+++ /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.h	(revision 26131)
@@ -22,30 +22,11 @@
 
 		int      nrows;
-		int      ncols;
-		bool     dofsymmetrical;
+		int      fsize;
+		int      ssize;
 		IssmDouble*  values;
 
-		//gset
 		int*     gglobaldoflist;
-
-		/*row wise: */
-		//fset
-		int      row_fsize;
-		int*     row_flocaldoflist;
-		int*     row_fglobaldoflist;
-		//sset
-		int      row_ssize;
-		int*     row_slocaldoflist;
-		int*     row_sglobaldoflist;
-
-		/*column wise: */
-		//fset
-		int      col_fsize;
-		int*     col_flocaldoflist;
-		int*     col_fglobaldoflist;
-		//sset
-		int      col_ssize;
-		int*     col_slocaldoflist;
-		int*     col_sglobaldoflist;
+		int*     fglobaldoflist;
+		int*     sglobaldoflist;
 
 		/*ElementMatrix constructors, destructors*/
Index: /issm/trunk-jpl/src/c/classes/matrix/ElementVector.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/matrix/ElementVector.cpp	(revision 26130)
+++ /issm/trunk-jpl/src/c/classes/matrix/ElementVector.cpp	(revision 26131)
@@ -21,7 +21,6 @@
 
 	this->nrows=0;
+	this->fsize=0;
 	this->values=NULL;
-	this->fsize=0;
-	this->flocaldoflist=NULL;
 	this->fglobaldoflist=NULL;
 
@@ -31,8 +30,6 @@
 
 	/*intermediaries*/
-	int i,j,counter;
+	int i,j;
 	int gsize,fsize;
-	int* P=NULL;
-	bool found;
 
 	/*If one of the two matrix is NULL, we copy the other one*/
@@ -50,10 +47,10 @@
 
 	/*Initialize itransformation matrix pe[P[i]] = pe2[i]*/
-	P=xNew<int>(pe2->nrows);
+	int* P=xNew<int>(pe2->nrows);
 
 	/*1: Get the new numbering of pe2 and get size of the new matrix*/
 	gsize=pe1->nrows;
 	for(i=0;i<pe2->nrows;i++){
-		found=false;
+		bool found=false;
 		for(j=0;j<pe1->nrows;j++){
 			if(pe2->gglobaldoflist[i]==pe1->gglobaldoflist[j]){
@@ -71,40 +68,20 @@
 	/*Gset and values*/
 	this->gglobaldoflist=xNew<int>(this->nrows);
+	this->fglobaldoflist=xNew<int>(this->nrows);
 	this->values=xNewZeroInit<IssmDouble>(this->nrows);
 	for(i=0;i<pe1->nrows;i++){
 		this->values[i] += pe1->values[i];
 		this->gglobaldoflist[i]=pe1->gglobaldoflist[i];
+		this->fglobaldoflist[i]=pe1->fglobaldoflist[i];
 	}
 	for(i=0;i<pe2->nrows;i++){
 		this->values[P[i]] += pe2->values[i];
 		this->gglobaldoflist[P[i]]=pe2->gglobaldoflist[i];
+		this->fglobaldoflist[P[i]]=pe2->fglobaldoflist[i];
 	}
 
 	/*Fset*/
-	fsize=pe1->fsize;
-	for(i=0;i<pe2->fsize;i++){
-		if(P[pe2->flocaldoflist[i]] >= pe1->nrows) fsize++;
-	}
-	this->fsize=fsize;
-	if(fsize){
-		this->flocaldoflist =xNew<int>(fsize);
-		this->fglobaldoflist=xNew<int>(fsize);
-		for(i=0;i<pe1->fsize;i++){
-			this->flocaldoflist[i] =pe1->flocaldoflist[i];
-			this->fglobaldoflist[i]=pe1->fglobaldoflist[i];
-		}
-		counter=pe1->fsize;
-		for(i=0;i<pe2->fsize;i++){
-			if(P[pe2->flocaldoflist[i]] >= pe1->nrows){
-				this->flocaldoflist[counter] =P[pe2->flocaldoflist[i]];
-				this->fglobaldoflist[counter]=pe2->fglobaldoflist[i];
-				counter++;
-			}
-		}
-	}
-	else{
-		this->flocaldoflist=NULL;
-		this->fglobaldoflist=NULL;
-	}
+	this->fsize=0;
+	for(i=0;i<this->nrows;i++) if(this->fglobaldoflist[i]>=0) this->fsize++;
 
 	/*clean-up*/
@@ -139,13 +116,10 @@
 	/*Get fsize*/
 	this->fsize=GetNumberOfDofs(nodes,numnodes,FsetEnum,approximation);
-	this->flocaldoflist =GetLocalDofList( nodes,numnodes,FsetEnum,approximation);
 	this->fglobaldoflist=GetGlobalDofList(nodes,numnodes,FsetEnum,approximation);
 }
 /*}}}*/
 ElementVector::~ElementVector(){/*{{{*/
-
 	xDelete<IssmDouble>(this->values);
 	xDelete<int>(this->gglobaldoflist);
-	xDelete<int>(this->flocaldoflist);
 	xDelete<int>(this->fglobaldoflist);
 }
@@ -155,21 +129,9 @@
 void ElementVector::AddToGlobal(Vector<IssmDouble>* pf){/*{{{*/
 
-	int i;
-	IssmDouble* localvalues=NULL;
-
 	/*In debugging mode, check consistency (no NaN, and values not too big)*/
 	this->CheckConsistency();
 
 	if(this->fsize){
-		/*first, retrieve values that are in the f-set from the g-set values vector: */
-		localvalues=xNew<IssmDouble>(this->fsize);
-		for(i=0;i<this->fsize;i++){
-			localvalues[i]=this->values[this->flocaldoflist[i]];
-		}
-		/*add local values into global  vector, using the fglobaldoflist: */
-		pf->SetValues(this->fsize,this->fglobaldoflist,localvalues,ADD_VAL);
-
-		/*Free ressources:*/
-		xDelete<IssmDouble>(localvalues);
+		pf->SetValues(this->nrows,this->fglobaldoflist,this->values,ADD_VAL);
 	}
 
@@ -193,8 +155,7 @@
 	_printf_("Element Vector echo:\n");
 	_printf_("   nrows: " << nrows << "\n");
+	_printf_("   fsize: " << fsize << "\n");
 	_printf_("   values:\n");
-	for(i=0;i<nrows;i++){
-		_printf_(setw(4) << right << i << ": " << setw(10) << values[i] << "\n");
-	}
+	for(i=0;i<nrows;i++) _printf_(setw(4) << right << i << ": " << setw(10) << values[i] << "\n");
 
 	_printf_("   gglobaldoflist (" << gglobaldoflist << "): ");
@@ -202,8 +163,4 @@
 	_printf_(" \n");
 
-	_printf_("   fsize: " << fsize << "\n");
-	_printf_("   flocaldoflist  (" << flocaldoflist << "): ");
-	if(flocaldoflist) for(i=0;i<fsize;i++) _printf_(" " << flocaldoflist[i] );
-	_printf_(" \n");
 	_printf_("   fglobaldoflist (" << fglobaldoflist << "): ");
 	if(fglobaldoflist) for(i=0;i<fsize;i++) _printf_(" " << fglobaldoflist[i] );
@@ -216,4 +173,5 @@
 
 	this->nrows =pe->nrows;
+	this->fsize =pe->fsize;
 
 	this->values=xNew<IssmDouble>(this->nrows);
@@ -223,30 +181,13 @@
 	xMemCpy<int>(this->gglobaldoflist,pe->gglobaldoflist,this->nrows);
 
-	this->fsize=pe->fsize;
+	this->fglobaldoflist=xNew<int>(this->nrows);
+	xMemCpy<int>(this->fglobaldoflist,pe->fglobaldoflist,this->nrows);
+}
+/*}}}*/
+void ElementVector::InsertIntoGlobal(Vector<IssmDouble>* pf){/*{{{*/
+
 	if(this->fsize){
-		this->flocaldoflist=xNew<int>(this->fsize);
-		xMemCpy<int>(this->flocaldoflist,pe->flocaldoflist,this->fsize);
-		this->fglobaldoflist=xNew<int>(this->fsize);
-		xMemCpy<int>(this->fglobaldoflist,pe->fglobaldoflist,this->fsize);
-	}
-	else{
-		this->flocaldoflist=NULL;
-		this->fglobaldoflist=NULL;
-	}
-}
-/*}}}*/
-void ElementVector::InsertIntoGlobal(Vector<IssmDouble>* pf){/*{{{*/
-
-	if(this->fsize){
-		/*first, retrieve values that are in the f-set from the g-set values vector: */
-		IssmDouble* localvalues=xNew<IssmDouble>(this->fsize);
-		for(int i=0;i<this->fsize;i++){
-			localvalues[i]=this->values[this->flocaldoflist[i]];
-		}
 		/*add local values into global  vector, using the fglobaldoflist: */
-		pf->SetValues(this->fsize,this->fglobaldoflist,localvalues,INS_VAL);
-
-		/*Free ressources:*/
-		xDelete<IssmDouble>(localvalues);
+		pf->SetValues(this->nrows,this->fglobaldoflist,this->values,INS_VAL);
 	}
 
@@ -280,5 +221,5 @@
 	/*Checks in debugging mode*/
 	_assert_(bsize>0 && bsize<this->nrows && this->values && Ke); 
-	_assert_(Ke->nrows==Ke->ncols && this->nrows==Ke->nrows);
+	_assert_(this->nrows==Ke->nrows);
 
 	/*Intermediaries*/
@@ -317,6 +258,6 @@
 	Fb  = xNew<IssmDouble>(bsize);
 	Fi  = xNew<IssmDouble>(isize);
-	for(i=0;i<isize;i++) for(j=0;j<bsize;j++) Kib[i*bsize+j] = Ke->values[iindices[i]*Ke->ncols + bindices[j]];
-	for(i=0;i<bsize;i++) for(j=0;j<bsize;j++) Kbb[i*bsize+j] = Ke->values[bindices[i]*Ke->ncols + bindices[j]];
+	for(i=0;i<isize;i++) for(j=0;j<bsize;j++) Kib[i*bsize+j] = Ke->values[iindices[i]*Ke->nrows + bindices[j]];
+	for(i=0;i<bsize;i++) for(j=0;j<bsize;j++) Kbb[i*bsize+j] = Ke->values[bindices[i]*Ke->nrows + bindices[j]];
 	for(i=0;i<bsize;i++) Fb[i] = this->values[bindices[i]];
 	for(i=0;i<isize;i++) Fi[i] = this->values[iindices[i]];
Index: /issm/trunk-jpl/src/c/classes/matrix/ElementVector.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/matrix/ElementVector.h	(revision 26130)
+++ /issm/trunk-jpl/src/c/classes/matrix/ElementVector.h	(revision 26131)
@@ -22,13 +22,8 @@
 	public:
 		int         nrows;
-		IssmDouble* values;
-
-		//gset
-		int* gglobaldoflist;
-
-		//fset
-		int  fsize;
-		int* flocaldoflist;
-		int* fglobaldoflist;
+		int         fsize;
+		IssmDouble *values;
+		int        *gglobaldoflist;
+		int        *fglobaldoflist;
 
 		/*ElementVector constructors, destructors*/
Index: /issm/trunk-jpl/src/c/toolkits/issm/Bucket.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/Bucket.h	(revision 26130)
+++ /issm/trunk-jpl/src/c/toolkits/issm/Bucket.h	(revision 26131)
@@ -34,37 +34,78 @@
 			this->Initialize();
 		} /*}}}*/
-		Bucket(int min,int* idxmin,int nin,int* idxnin,doubletype* valuesin,InsMode modein){ /*{{{*/
+		Bucket(int m_in,int* idxm_in,int n_in,int* idxn_in,doubletype* values_in,InsMode mode_in){ /*{{{*/
 
 			this->Initialize();
 
+			/*Check whether we have negative indices?*/
+			int m_numneg = 0;
+			int n_numneg = 0;
+			for(int i=0;i<m_in;i++) if(idxm_in[i]<0) m_numneg++;
+			for(int i=0;i<n_in;i++) if(idxn_in[i]<0) n_numneg++;
+
 			this->type=MATRIX_BUCKET;
-			this->m=min;
-			this->n=nin;
-			this->mode=modein;
+			this->m=m_in-m_numneg;
+			this->n=n_in-n_numneg;
+			this->mode=mode_in;
+			if(this->m*this->n){
+				this->idxm=xNew<int>(this->m); 
+				this->idxn=xNew<int>(this->n); 
+				this->values=xNew<doubletype>(this->n*this->m);
+
+				if(m_numneg==0 && n_numneg==0){
+					xMemCpy(this->values,values_in,this->n*this->m);
+					xMemCpy(this->idxm,idxm_in,this->m);
+					xMemCpy(this->idxn,idxn_in,this->n);
+				}
+				else{
+					int count;
+					count =0; for(int i=0;i<m_in;i++) if(idxm_in[i]>=0) this->idxm[count++] = idxm_in[i];
+					count =0; for(int i=0;i<n_in;i++) if(idxn_in[i]>=0) this->idxn[count++] = idxn_in[i];
+					count = 0;
+					for(int i=0;i<m_in;i++){
+						if(idxm_in[i]>=0){
+							for(int j=0;j<n_in;j++){
+								if(idxn_in[j]>=0){
+									this->values[count++] = values_in[i*n_in+j];
+								}
+							}
+						}
+					}
+				}
+			}
+			else{
+				/*Set both as 0 to save time*/
+				this->m = 0;
+				this->n = 0;
+			}
+		} /*}}}*/
+		Bucket(int m_in,int* idxm_in,doubletype* values_in,InsMode mode_in){ /*{{{*/ 
+			this->Initialize();
+
+			/*Check whether we have negative indices?*/
+			int numneg = 0;
+			for(int i=0;i<m_in;i++) if(idxm_in[i]<0) numneg++;
+
+			this->type=VECTOR_BUCKET; 
+			this->m=m_in-numneg;
+			this->mode=mode_in;
 			if(this->m){
-				this->idxm=xNew<int>(this->m); 
-				xMemCpy(this->idxm,idxmin,this->m);
-			}
-			if(this->n){
-				this->idxn=xNew<int>(this->n); 
-				xMemCpy(this->idxn,idxnin,this->n);
-			}
-			if(this->m*this->n){
-				this->values=xNew<doubletype>(this->n*this->m);
-				xMemCpy(this->values,valuesin,this->n*this->m);
-			}
-		} /*}}}*/
-		Bucket(int min,int* idxmin,doubletype* valuesin,InsMode modein){ /*{{{*/ 
-			this->Initialize();
-
-			this->type=VECTOR_BUCKET; 
-			this->m=min;
-			this->mode=modein;
-			if(this->m){
-				this->idxm=xNew<int>(this->m); 
-				xMemCpy(this->idxm,idxmin,this->m);
-
+				this->idxm=xNew<int>(this->m);
 				this->values=xNew<doubletype>(this->m);
-				xMemCpy(this->values,valuesin,this->m);
+
+				if(numneg==0){
+					xMemCpy(this->idxm,idxm_in,this->m);
+					xMemCpy(this->values,values_in,this->m);
+				}
+				else{
+					int count = 0;
+					for(int i=0;i<m_in;i++){
+						if(idxm_in[i]>=0){
+							this->idxm[count]   = idxm_in[i];
+							this->values[count] = values_in[i];
+							count++;
+						}
+					}
+				}
 			}
 		} /*}}}*/
@@ -88,5 +129,5 @@
 		/*object virtual functions definitions:*/
 		void    Echo(){ /*{{{*/
-			_printf_("Bucket echo (cpu #: "<<IssmComm::GetRank()<<")\n");
+			_printf_("Bucket echo (cpu #"<<IssmComm::GetRank()<<")\n");
 			_printf_("bucket type: " << type << "\n");
 			_printf_("num rows: "<<this->m<<" num cols: "<<this->n << "\n");
@@ -95,5 +136,5 @@
 			int i,j;
 
-			_printf_("Bucket echo (cpu #: "<<IssmComm::GetRank()<<")\n");
+			_printf_("Bucket echo (cpu #"<<IssmComm::GetRank()<<")\n");
 			_printf_("bucket type: " << type << "\n");
 			_printf_("num rows: "<<this->m<<" num cols: "<<this->n << "\n");
Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp	(revision 26130)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp	(revision 26131)
@@ -37,4 +37,5 @@
 	VecSetSizes(this->vector,m,M);
 	VecSetFromOptions(this->vector);
+	VecSetOption(this->vector,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
 	VecSet(this->vector,0.0);
 
Index: /issm/trunk-jpl/src/c/toolkits/petsc/patches/NewMat.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/patches/NewMat.cpp	(revision 26130)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/patches/NewMat.cpp	(revision 26131)
@@ -111,4 +111,5 @@
 	MatSetSizes(outmatrix,m,n,M,N);
 	MatSetFromOptions(outmatrix);
+	MatSetOption(outmatrix,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
 
 	/*preallocation  according to type: */
Index: /issm/trunk-jpl/src/c/toolkits/petsc/patches/NewVec.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/patches/NewVec.cpp	(revision 26130)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/patches/NewVec.cpp	(revision 26131)
@@ -35,4 +35,5 @@
 	VecSetSizes(vector,local_size,PETSC_DECIDE); 
 	VecSetFromOptions(vector); 
+	VecSetOption(vector,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
 
 	return vector;