Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 19106)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 19107)
@@ -123,9 +123,327 @@
 }/*}}}*/
 ElementMatrix* LevelsetAnalysis::CreateKMatrix(Element* element){/*{{{*/
-	_error_("Not implemented yet");
+
+	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;
 }/*}}}*/
 ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/
 	
-	_error_("Not implemented yet");
+	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;
 }/*}}}*/
 void           LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/