Index: /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 24432)
+++ /issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 24433)
@@ -82,6 +82,4 @@
 }/*}}}*/
 ElementMatrix* ExtrapolationAnalysis::CreateKMatrix(Element* element){/*{{{*/
-
-	if(!element->AnyFSet()) return NULL;
 
 	/*Intermediaries */
Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 24432)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 24433)
@@ -789,61 +789,4 @@
 	return pe;
 }/*}}}*/
-void           MasstransportAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
-	 * For node i, Bi can be expressed in the actual coordinate system
-	 * by:
-	 *       Bi=[ N ]
-	 *          [ N ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions*/
-	IssmDouble* basis=xNew<IssmDouble>(numnodes);
-	element->NodalFunctions(basis,gauss);
-
-	/*Build B: */
-	for(int i=0;i<numnodes;i++){
-		for(int j=0;j<dim;j++){
-			B[numnodes*j+i] = basis[i];
-		}
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
-}/*}}}*/
-void           MasstransportAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
-	 * For node i, Bi' can be expressed in the actual coordinate system
-	 * by:
-	 *       Bi_prime=[ dN/dx ]
-	 *                [ dN/dy ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
-	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B': */
-	for(int i=0;i<numnodes;i++){
-		for(int j=0;j<dim;j++){
-			Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
-		}
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-
-}/*}}}*/
 void           MasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h	(revision 24432)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h	(revision 24433)
@@ -30,6 +30,4 @@
 		ElementVector* CreatePVectorCG(Element* element);
 		ElementVector* CreatePVectorDG(Element* element);
-		void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void           GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
 		void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
 		void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 24432)
+++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 24433)
@@ -43,5 +43,5 @@
 		for (i=0;i<femmodel->elements->Size();i++){
 			element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-			//if(!element->AnyFSet()) continue;
+			if(!element->AnyFSet() && analysisenum!=StressbalanceAnalysisEnum) continue;
 			ElementMatrix* Ke = analysis->CreateKMatrix(element);
 			ElementVector* pe = analysis->CreatePVector(element);
@@ -77,5 +77,5 @@
 	for (i=0;i<femmodel->elements->Size();i++){
 		element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-		//if(!element->AnyFSet()) continue;
+		if(!element->AnyFSet() && analysisenum!=StressbalanceAnalysisEnum) continue;
 		ElementMatrix* Ke = analysis->CreateKMatrix(element);
 		ElementVector* pe = analysis->CreatePVector(element);