Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 17800)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 17801)
@@ -89,170 +89,6 @@
 ElementMatrix* LevelsetAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
-	if(!element->IsOnBase()) return NULL;
-	Element* basalelement = element->SpawnBasalElement();
-
-	/*Intermediaries */
-	int  dim, domaintype;
-	int i, row, col;
-	IssmDouble kappa;
-	IssmDouble Jdet, dt, D_scalar;
-	IssmDouble h,hx,hy,hz;
-	IssmDouble vel, calvingrate;
-	IssmDouble norm_dlsf;
-	IssmDouble* xyz_list = NULL;
-
-	/*Get problem dimension*/
-	basalelement->FindParam(&domaintype,DomainTypeEnum);
-	switch(domaintype){
-		case Domain2DverticalEnum:   dim = 1; break;
-		case Domain2DhorizontalEnum: dim = 2; break;
-		case Domain3DEnum:           dim = 2; break;
-		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-	}
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes    = basalelement->GetNumberOfNodes();
-
-	/*Initialize Element vector and other vectors*/
-	ElementMatrix* Ke       = basalelement->NewElementMatrix();
-	IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
-	IssmDouble*    B        = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble*    Bprime   = xNew<IssmDouble>(dim*numnodes);
-	IssmDouble*    D        = xNew<IssmDouble>(dim*dim);
-	IssmDouble*    v        = xNew<IssmDouble>(dim);
-	IssmDouble*    w        = xNew<IssmDouble>(dim);
-	IssmDouble*    c        = xNew<IssmDouble>(dim);
-	IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
-
-	/*Retrieve all inputs and parameters*/
-	basalelement->GetVerticesCoordinates(&xyz_list);
-	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
-	Input* vx_input=NULL;
-	Input* vy_input=NULL;
-	if(domaintype==Domain2DhorizontalEnum){
-		vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
-		vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
-	}
-	else{
-		if(dim==1){
-			vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
-		}
-		if(dim==2){
-			vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
-			vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
-		}
-	}
-
-	Input* lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
-	Input* lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
-	Input* calvingrate_input  = basalelement->GetInput(MasstransportCalvingrateEnum);     _assert_(calvingrate_input);
-	
-	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-
-		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		D_scalar=gauss->weight*Jdet;
-
-		/* Transient */
-		if(dt!=0.){
-			basalelement->NodalFunctions(basis,gauss);
-			TripleMultiply(basis,numnodes,1,0,
-						&D_scalar,1,1,0,
-						basis,1,numnodes,0,
-						&Ke->values[0],1);
-			D_scalar*=dt;
-		}
-
-		/* Advection */
-		GetB(B,basalelement,xyz_list,gauss); 
-		GetBprime(Bprime,basalelement,xyz_list,gauss); 
-		vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity 
-		vy_input->GetInputValue(&v[1],gauss); 
-		lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
-		lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
-		calvingrate_input->GetInputValue(&calvingrate,gauss);
-
-		norm_dlsf=0.;
-		for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
-		norm_dlsf=sqrt(norm_dlsf);
-
-		if(norm_dlsf>1.e-10)
-			for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
-		else
-			for(i=0;i<dim;i++) c[i]=0.;
-		
-		for(i=0;i<dim;i++) w[i]=v[i]-c[i];
-
-		for(row=0;row<dim;row++)
-			for(col=0;col<dim;col++)
-				if(row==col)
-					D[row*dim+col]=D_scalar*w[row];
-				else
-				   D[row*dim+col]=0.;
-
-		TripleMultiply(B,dim,numnodes,1,
-					D,dim,dim,0,
-					Bprime,dim,numnodes,0,
-					&Ke->values[0],1);
-
-		/* Stabilization */
-		int stabilization=2;
-		vel=0.;
-		for(i=0;i<dim;i++) vel+=v[i]*v[i];
-		vel=sqrt(vel)+1.e-14;
-		switch(stabilization){
-			case 0:
-				// no stabilization, do nothing
-				break;
-			case 1:
-				/* Artificial Diffusion */
-				basalelement->ElementSizes(&hx,&hy,&hz);
-				h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) ); 
-				kappa=h*vel/2.;
-				for(row=0;row<dim;row++)
-					for(col=0;col<dim;col++)
-					if(row==col)
-						D[row*dim+col]=D_scalar*kappa;
-					else
-						D[row*dim+col]=0.;
-
-				TripleMultiply(Bprime,dim,numnodes,1,
-							D,dim,dim,0,
-							Bprime,dim,numnodes,0,
-							&Ke->values[0],1);
-				break;	
-			case 2:
-				/* Streamline Upwinding */
-				basalelement->ElementSizes(&hx,&hy,&hz);
-				h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) );
-				for(row=0;row<dim;row++) 
-					for(col=0;col<dim;col++) 
-						D[row*dim+col] = D_scalar*h/(2.*vel)*v[row]*v[col];
-
-				TripleMultiply(Bprime,dim,numnodes,1,
-							D,dim,dim,0,
-							Bprime,dim,numnodes,0,
-							&Ke->values[0],1);
-				break;
-			default:
-				_error_("unknown type of stabilization in LevelsetAnalysis.cpp");
-		}
-	}
-
-	/*Clean up and return*/
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(basis);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(D);
-	xDelete<IssmDouble>(Bprime);
-	xDelete<IssmDouble>(v);
-	xDelete<IssmDouble>(w);
-	xDelete<IssmDouble>(c);
-	xDelete<IssmDouble>(dlsf);
-	delete gauss;
-	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
-	return Ke;
+	_error_("NOT IMPLEMENTED YET");
+
 }/*}}}*/
 ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/