Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 19100)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 19101)
@@ -123,327 +123,9 @@
 }/*}}}*/
 ElementMatrix* LevelsetAnalysis::CreateKMatrix(Element* element){/*{{{*/
-
-	if(!element->IsOnBase()) return NULL;
-	Element* basalelement = element->SpawnBasalElement();
-
-	/*Intermediaries */
-	int  stabilization,dim, domaintype, calvinglaw;
-	bool iscalving;
-	int i, row, col;
-	IssmDouble kappa;
-	IssmDouble Jdet, dt, D_scalar;
-	IssmDouble h,hx,hy,hz;
-	IssmDouble vel;
-	IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate;
-	IssmDouble* xyz_list = NULL;
-
-	/*Get problem dimension and whether there is calving or not*/
-	basalelement->FindParam(&iscalving,TransientIscalvingEnum);
-	basalelement->FindParam(&domaintype,DomainTypeEnum);
-	basalelement->FindParam(&calvinglaw,CalvingLawEnum);
-	basalelement->FindParam(&stabilization,LevelsetStabilizationEnum);
-	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        = xNewZeroInit<IssmDouble>(dim);
-	IssmDouble*    m        = xNewZeroInit<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;
-	Input* calvingratex_input = NULL;
-	Input* calvingratey_input = NULL;
-	Input* lsf_slopex_input   = NULL;
-	Input* lsf_slopey_input   = NULL;
-	Input* calvingrate_input  = NULL;
-	Input* meltingrate_input  = NULL;
-
-	/*Load velocities*/
-	switch(domaintype){
-		case Domain2DverticalEnum:
-			vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
-			break;
-		case Domain2DhorizontalEnum:
-			vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
-			vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
-			break;
-		case Domain3DEnum:
-			vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
-			vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
-			break;
-		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-	}
-
-	/*Load calving inputs*/
-	if(iscalving){
-		switch(calvinglaw){
-			case DefaultCalvingEnum:
-				lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
-				if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
-				calvingrate_input = basalelement->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
-				meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
-				break;
-			case CalvingLevermannEnum:
-				switch(domaintype){
-					case Domain2DverticalEnum:
-						calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
-						break;
-					case Domain2DhorizontalEnum:
-						calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
-						calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
-						break;
-					case Domain3DEnum:
-						calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
-						calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
-						break;
-					default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-				}
-				meltingrate_input = basalelement->GetInput(CalvinglevermannMeltingrateEnum);     _assert_(meltingrate_input);
-				break;
-			case CalvingPiEnum:
-				switch(domaintype){
-					case Domain2DverticalEnum:
-						calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
-						break;
-					case Domain2DhorizontalEnum:
-						calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
-						calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
-						break;
-					case Domain3DEnum:
-						calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
-						calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
-						break;
-					default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-				}
-				meltingrate_input = basalelement->GetInput(CalvingpiMeltingrateEnum);     _assert_(meltingrate_input);
-				break;
-			case CalvingDevEnum:
-				switch(domaintype){
-					case Domain2DverticalEnum:
-						calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
-						break;
-					case Domain2DhorizontalEnum:
-						calvingratex_input=basalelement->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
-						calvingratey_input=basalelement->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
-						break;
-					case Domain3DEnum:
-						calvingratex_input=basalelement->GetInput(CalvingratexAverageEnum); _assert_(calvingratex_input);
-						calvingratey_input=basalelement->GetInput(CalvingrateyAverageEnum); _assert_(calvingratey_input);
-						break;
-					default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-				}
-				meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
-				break;
-			default:
-				_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
-		}
-	}
-
-	/* 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);
-		vy_input->GetInputValue(&v[1],gauss); 
-
-		/*Get calving speed*/
-		if(iscalving){
-			switch(calvinglaw){
-				case DefaultCalvingEnum:
-					lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
-					if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
-					calvingrate_input->GetInputValue(&calvingrate,gauss);
-					meltingrate_input->GetInputValue(&meltingrate,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; m[i]=meltingrate*dlsf[i]/norm_dlsf;
-					 }
-					else
-					 for(i=0;i<dim;i++){
-						 c[i]=0.; m[i]=0.;
-					 }
-					break;
-
-				case CalvingLevermannEnum:
-				case CalvingPiEnum:
-				case CalvingDevEnum:
-					calvingratex_input->GetInputValue(&c[0],gauss);
-					if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
-					meltingrate_input->GetInputValue(&meltingrate,gauss);
-					norm_calving=0.;
-					for(i=0;i<dim;i++) norm_calving+=pow(c[i],2);
-					norm_calving=sqrt(norm_calving)+1.e-14;
-					for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
-					break;
-
-				default:
-					_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
-			}
-		}
-
-		/*Levelset speed is ice velocity - calving rate*/
-		for(i=0;i<dim;i++) w[i]=v[i]-c[i]-m[i];
-
-		/*Compute D*/
-		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 */
-		vel=0.;
-		for(i=0;i<dim;i++) vel+=w[i]*w[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*w[0]/vel,2) + pow(hy*w[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*w[0]/vel,2) + pow(hy*w[1]/vel,2) );
-				for(row=0;row<dim;row++) 
-					for(col=0;col<dim;col++) 
-						D[row*dim+col] = D_scalar*h/(2.*vel)*w[row]*w[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>(m);
-	xDelete<IssmDouble>(dlsf);
-	delete gauss;
-	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
-	return Ke;
+	_error_("Not implemented yet");
 }/*}}}*/
 ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/
 	
-	if(!element->IsOnBase()) return NULL;
-	Element* basalelement = element->SpawnBasalElement();
-
-	/*Intermediaries */
-	int i, ig, domaintype;
-	IssmDouble  Jdet,dt;
-	IssmDouble  lsf;
-	IssmDouble* xyz_list = NULL;
-	
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = basalelement->GetNumberOfNodes();
-
-	/*Initialize Element vector*/
-	ElementVector* pe = basalelement->NewElementVector();
-	basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
-	
-	if(dt!=0.){
-		/*Initialize basis vector*/
-		IssmDouble*    basis = xNew<IssmDouble>(numnodes);
-
-		/*Retrieve all inputs and parameters*/
-		basalelement->GetVerticesCoordinates(&xyz_list);
-		Input* levelset_input     = basalelement->GetInput(MaskIceLevelsetEnum);                    _assert_(levelset_input);
-
-		/* Start  looping on the number of gaussian points: */
-		Gauss* gauss=basalelement->NewGauss(2);
-		for(ig=gauss->begin();ig<gauss->end();ig++){
-			gauss->GaussPoint(ig);
-
-			basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
-			basalelement->NodalFunctions(basis,gauss);
-
-			/* old function value */
-			levelset_input->GetInputValue(&lsf,gauss);
-			for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*lsf*basis[i];
-		}
-
-		/*Clean up and return*/
-		xDelete<IssmDouble>(xyz_list);
-		xDelete<IssmDouble>(basis);
-		basalelement->FindParam(&domaintype,DomainTypeEnum);
-		if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
-		delete gauss;
-	}
-
-	return pe;
+	_error_("Not implemented yet");
 }/*}}}*/
 void           LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/