Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp	(revision 20664)
+++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp	(revision 20665)
@@ -49,108 +49,4 @@
 ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
-	/*compute all stiffness matrices for this element*/
-	ElementMatrix* Ke1=CreateKMatrixVolume(element);
-	ElementMatrix* Ke2=CreateKMatrixSurface(element);
-	ElementMatrix* Ke3=CreateKMatrixBed(element);
-	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
-
-	/*clean-up and return*/
-	delete Ke1;
-	delete Ke2;
-	delete Ke3;
-	return Ke;
-}/*}}}*/
-ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixBed(Element* element){/*{{{*/
-
-	if(!element->IsOnBase()) return NULL;
-
-	/*Intermediaries */
-	int         dim;
-	IssmDouble  Jdet,D,normal[3];
-	IssmDouble *xyz_list_base = NULL;
-
-	/*Get dimension*/
-	element->FindParam(&dim,DomainDimensionEnum);
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Initialize Element vector and other vectors*/
-	ElementMatrix* Ke    = element->NewElementMatrix();
-	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
-
-	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinatesBase(&xyz_list_base);
-	element->NormalBase(&normal[0],xyz_list_base);
-
-	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGaussBase(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-
-		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-		element->NodalFunctions(basis,gauss);
-		D = - gauss->weight*Jdet*normal[dim-1]; 
-
-		TripleMultiply(basis,1,numnodes,1,
-					&D,1,1,0,
-					basis,1,numnodes,0,
-					&Ke->values[0],1);
-	}
-
-	/*Clean up and return*/
-	xDelete<IssmDouble>(xyz_list_base);
-	xDelete<IssmDouble>(basis);
-	delete gauss;
-	return Ke;
-}
-/*}}}*/
-ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/
-
-	if(!element->IsOnSurface()) return NULL;
-
-	/*Intermediaries */
-	int         dim;
-	IssmDouble  Jdet,D,normal[3];
-	IssmDouble *xyz_list_top = NULL;
-
-	/*Get dimension*/
-	element->FindParam(&dim,DomainDimensionEnum);
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Initialize Element vector and other vectors*/
-	ElementMatrix* Ke    = element->NewElementMatrix();
-	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
-
-	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinatesTop(&xyz_list_top);
-	element->NormalTop(&normal[0],xyz_list_top);
-
-	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGaussTop(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-
-		element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
-		element->NodalFunctions(basis,gauss);
-		D = - gauss->weight*Jdet*normal[dim-1]; 
-
-		TripleMultiply(basis,1,numnodes,1,
-					&D,1,1,0,
-					basis,1,numnodes,0,
-					&Ke->values[0],1);
-	}
-
-	/*Clean up and return*/
-	xDelete<IssmDouble>(xyz_list_top);
-	xDelete<IssmDouble>(basis);
-	delete gauss;
-	return Ke;
-}
-/*}}}*/
-ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
-
 	/*Intermediaries */
 	IssmDouble  Jdet,D;
@@ -166,6 +62,5 @@
 	/*Initialize Element vector and other vectors*/
 	ElementMatrix* Ke     = element->NewElementMatrix();
-	IssmDouble*    B      = xNew<IssmDouble>(numnodes);
-	IssmDouble*    Bprime = xNew<IssmDouble>(numnodes);
+	IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -178,45 +73,22 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		element->NodalFunctions(Bprime,gauss);
-		GetB(B,element,dim,xyz_list,gauss);
 		D=gauss->weight*Jdet;
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
-		TripleMultiply(B,1,numnodes,1,
-					&D,1,1,0,
-					Bprime,1,numnodes,0,
-					&Ke->values[0],1);
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += gauss->weight*Jdet*(dbasis[(dim-1)*numnodes+i]*dbasis[(dim-1)*numnodes+j]);
+			}
+		}
 	} 
 
 	/*Clean up and return*/
 	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(Bprime);
 	delete gauss;
 	return Ke;
-}
-/*}}}*/
+}/*}}}*/
 ElementVector* ExtrudeFromBaseAnalysis::CreatePVector(Element* element){/*{{{*/
 	return NULL;
 }/*}}}*/
-void           ExtrudeFromBaseAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*	Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
-		where hi is the interpolation function for node i.*/
-
-	/*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++){
-		B[i] = dbasis[(dim-1)*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-}
-/*}}}*/
 void           ExtrudeFromBaseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	   _error_("not implemented yet");
Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h	(revision 20664)
+++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h	(revision 20665)
@@ -25,9 +25,5 @@
 		ElementMatrix* CreateJacobianMatrix(Element* element);
 		ElementMatrix* CreateKMatrix(Element* element);
-		ElementMatrix* CreateKMatrixBed(Element* element);
-		ElementMatrix* CreateKMatrixSurface(Element* element);
-		ElementMatrix* CreateKMatrixVolume(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void           GetB(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/classes/Materials/Matice.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp	(revision 20664)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp	(revision 20665)
@@ -642,4 +642,5 @@
 	/*Get viscosity*/
 	this->GetViscosity(&viscosity,eps_eff);
+	_assert_(!xIsNan<IssmDouble>(viscosity));
 
 	/*Assign output pointer*/
Index: /issm/trunk-jpl/src/c/classes/Vertex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Vertex.cpp	(revision 20664)
+++ /issm/trunk-jpl/src/c/classes/Vertex.cpp	(revision 20665)
@@ -164,4 +164,5 @@
 			this->y = newy;
 			vy->SetValue(this->pid,vely,INS_VAL);
+			_assert_(!xIsNan<IssmDouble>(vely));
 			return;
 		case Domain3DEnum:
@@ -171,4 +172,5 @@
 			this->z = newz;
 			vz->SetValue(this->pid,velz,INS_VAL);
+			_assert_(!xIsNan<IssmDouble>(velz));
 			return;
 		default: