Index: /issm/trunk-jpl/src/c/CMakeLists.txt
===================================================================
--- /issm/trunk-jpl/src/c/CMakeLists.txt	(revision 16230)
+++ /issm/trunk-jpl/src/c/CMakeLists.txt	(revision 16231)
@@ -515,6 +515,4 @@
 					./toolkits/petsc/patches/NewMat.cpp
 					./toolkits/petsc/patches/VecFree.cpp
-					./toolkits/petsc/patches/PetscMatrixToDoubleMatrix.cpp
-					./toolkits/petsc/patches/PetscVectorToDoubleVector.cpp
 					./toolkits/petsc/patches/KSPFree.cpp
 					./toolkits/petsc/patches/MatFree.cpp
Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 16230)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 16231)
@@ -767,6 +767,4 @@
 					./toolkits/petsc/patches/NewMat.cpp\
 					./toolkits/petsc/patches/VecFree.cpp\
-					./toolkits/petsc/patches/PetscMatrixToDoubleMatrix.cpp\
-					./toolkits/petsc/patches/PetscVectorToDoubleVector.cpp\
 					./toolkits/petsc/patches/KSPFree.cpp\
 					./toolkits/petsc/patches/MatFree.cpp\
Index: /issm/trunk-jpl/src/c/bamg/Curve.cpp
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Curve.cpp	(revision 16230)
+++ /issm/trunk-jpl/src/c/bamg/Curve.cpp	(revision 16231)
@@ -20,11 +20,4 @@
 
 	/*Methods*/
-	/*FUNCTION Curve::Reverse {{{*/
-	void Curve::Reverse() {
-		/*reverse the direction of the curve */
-		Exchange(FirstEdge,LastEdge);
-		Exchange(FirstVertexIndex,LastVertexIndex);
-	}
-	/*}}}*/
 	/*FUNCTION Curve::Set {{{*/
 	void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew){
Index: /issm/trunk-jpl/src/c/bamg/Curve.h
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Curve.h	(revision 16230)
+++ /issm/trunk-jpl/src/c/bamg/Curve.h	(revision 16231)
@@ -19,5 +19,4 @@
 			//Methods
 			Curve();
-			void Reverse(void);
 			void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew);
 	};
Index: /issm/trunk-jpl/src/c/bamg/GeomEdge.cpp
===================================================================
--- /issm/trunk-jpl/src/c/bamg/GeomEdge.cpp	(revision 16230)
+++ /issm/trunk-jpl/src/c/bamg/GeomEdge.cpp	(revision 16231)
@@ -175,8 +175,4 @@
 		type |= 64;/*=>6th digit to 1*/ 
 	}/*}}}*/
-	  /*FUNCTION GeomEdge::Tg{{{*/
-	int    GeomEdge::Tg(int i) const  {
-		return i==0 ? TgA() : TgB();
-	}/*}}}*/
 	/*FUNCTION GeomEdge::TgA{{{*/
 	int    GeomEdge::TgA()     const  {
Index: /issm/trunk-jpl/src/c/bamg/GeomEdge.h
===================================================================
--- /issm/trunk-jpl/src/c/bamg/GeomEdge.h	(revision 16230)
+++ /issm/trunk-jpl/src/c/bamg/GeomEdge.h	(revision 16231)
@@ -28,5 +28,4 @@
 			R2     F(double theta) const ; // parametrization of the curve edge
 			double R1tg(double theta,R2 &t) const ; // 1/radius of curvature + tangente
-			int    Tg(int i) const;
 			int    Cracked() const;
 			int    TgA()     const;
Index: /issm/trunk-jpl/src/c/bamg/Mesh.cpp
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Mesh.cpp	(revision 16230)
+++ /issm/trunk-jpl/src/c/bamg/Mesh.cpp	(revision 16231)
@@ -3012,9 +3012,4 @@
 	}
 	/*}}}*/
-	/*FUNCTION Mesh::isCracked{{{*/
-	int Mesh::isCracked() const {
-		return NbCrackedVertices != 0;
-	}
-	/*}}}*/
 	/*FUNCTION Mesh::MakeGeomEdgeToEdge{{{*/
 	Edge** Mesh::MakeGeomEdgeToEdge() {
@@ -3817,78 +3812,4 @@
 	}
 	/*}}}*/
-	/*FUNCTION Mesh::VerticesRenumber{{{*/
-	void Mesh::VerticesRenumber(long * renu) {
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingVertex)*/
-
-		// warning be carfull because pointer
-		// from on mesh to over mesh 
-		//  --  so do ReNumbering at the beginning
-		BamgVertex * ve = vertices+nbv;
-		long it,ie,i;
-
-		_printf_("renumbering triangles\n");
-		for ( it=0;it<nbt;it++) 
-		 triangles[it].Renumbering(vertices,ve,renu);
-
-		_printf_("renumbering edges\n");
-		for ( ie=0;ie<nbe;ie++) 
-		 edges[ie].Renumbering(vertices,ve,renu);
-
-		_printf_("renumbering vertices on geom\n");
-		for (i=0;i< NbVerticesOnGeomVertex;i++)
-		  {
-			BamgVertex *v = VerticesOnGeomVertex[i].meshvertex;
-			if (v>=vertices && v < ve)
-			 VerticesOnGeomVertex[i].meshvertex=vertices+renu[GetId(v)];
-		  }
-
-		_printf_("renumbering vertices on edge\n");
-		for (i=0;i< NbVerticesOnGeomEdge;i++)
-		  {
-			BamgVertex *v =VerticesOnGeomEdge[i].meshvertex;
-			if (v>=vertices && v < ve)
-			 VerticesOnGeomEdge[i].meshvertex=vertices+renu[GetId(v)];
-		  }
-
-		_printf_("renumbering vertices on Bth vertex\n");
-		for (i=0;i< NbVertexOnBThVertex;i++)
-		  {
-			BamgVertex *v=VertexOnBThVertex[i].v;
-			if (v>=vertices && v < ve)
-			 VertexOnBThVertex[i].v=vertices+renu[GetId(v)];
-		  }
-
-		for (i=0;i< NbVertexOnBThEdge;i++)
-		  {
-			BamgVertex *v=VertexOnBThEdge[i].v;
-			if (v>=vertices && v < ve)
-			 VertexOnBThEdge[i].v=vertices+renu[GetId(v)];
-		  }
-
-		// move the Vertices without a copy of the array 
-		// be carefull not trivial code 
-		long j;
-		for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu
-		 if (renu[it] >= 0) // a new sub cycle
-			{ 
-			 i=it;
-			 BamgVertex ti=vertices[i],tj;
-			 while ( (j=renu[i]) >= 0){
-				 // i is old, and j is new 
-				 renu[i] = -1-renu[i]; // mark 
-				 tj = vertices[j];     // save new
-				 vertices[j]= ti;      // new <- old
-				 i=j;     // next 
-				 ti = tj;
-				}  
-			}
-		if (quadtree){
-			delete quadtree;
-			quadtree = new BamgQuadtree(this);
-		}
-
-		for ( it=0;it<nbv;it++) renu[i]= -renu[i]-1;
-	}
-	/*}}}*/
 /*FUNCTION Mesh::SetIntCoor{{{*/
 void Mesh::SetIntCoor(const char * strfrom) {
@@ -3951,113 +3872,4 @@
 
 	if (number_of_errors) _error_("Fatal error: some triangles have negative areas, see above");
-}
-/*}}}*/
-/*FUNCTION Mesh::ShowRegulaty{{{*/
-void  Mesh::ShowRegulaty() const {
-	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr>*/
-
-	const  double  sqrt32=sqrt(3.)*0.5; 
-	const double  aireKh=sqrt32*0.5;
-	D2  Beq(1,0),Heq(0.5,sqrt32);
-	D2xD2 Br(D2xD2(Beq,Heq).t());
-	D2xD2 B1r(Br.inv());
-	double gammamn=1e100,hmin=1e100;
-	double gammamx=0,hmax=0;
-	double beta=1e100;
-	double beta0=0;
-	double  alpha2=0;
-	double area=0,Marea=0;
-	// double cf= double(coefIcoor);
-	// double cf2= 6.*cf*cf;
-	int nt=0;
-	for (int it=0;it<nbt;it++)
-	 if ( triangles[it].link){
-		 Triangle &K=triangles[it];
-		 double  area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.;
-		 area+= area3;
-		 D2xD2 B_Kt(K[0],K[1],K[2]);
-		 D2xD2 B_K(B_Kt.t());
-		 D2xD2 B1K = Br*B_K.inv();
-		 D2xD2 BK =  B_K*B1r;
-		 D2xD2 B1B1 = B1K.t()*B1K;
-		 Metric MK(B1B1.x.x,B1B1.x.y,B1B1.y.y);
-		 EigenMetric VMK(MK);
-		 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1));
-		 double betaK=0;
-
-		 for (int j=0;j<3;j++)
-			{
-			 double he= Norme2(R2(K[j],K[(j+1)%3]));
-			 hmin=Min(hmin,he);
-			 hmax=Max(hmax,he);
-			 BamgVertex & v=K[j];
-			 D2xD2 M((Metric)v);
-			 betaK += sqrt(M.det());
-
-			 D2xD2 BMB = BK.t()*M*BK;
-			 Metric M1(BMB.x.x,BMB.x.y,BMB.y.y);
-			 EigenMetric VM1(M1);
-			 gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2);
-			 gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);		
-			}
-		 betaK *= area3;//  1/2 (somme sqrt(det))* area(K)
-		 Marea+= betaK;
-		 beta=min(beta,betaK);
-		 beta0=max(beta0,betaK);
-		}   
-	area*=3; 
-	gammamn=sqrt(gammamn);
-	gammamx=sqrt(gammamx);    
-	_printf_("   Adaptmesh info:\n");
-	_printf_("      number of triangles = " << nt << "\n");
-	_printf_("      hmin = " << hmin << ", hmax=" << hmax << "\n");
-	_printf_("      area = " << area << ", M area = " << Marea << ", M area/( |Khat| nt) = " <<  Marea/(aireKh*nt) << "\n");
-	_printf_("      infinite-regularity(?): min = " << gammamn << ", max = " << gammamx << "\n");
-	_printf_("      anisomax = " << pow(alpha2,0.5) << ", beta max = " << 1./pow(beta/aireKh,0.5) << ", min = " << 1./pow(beta0/aireKh,0.5) << "\n");
-}
-/*}}}*/
-/*FUNCTION Mesh::ShowHistogram{{{*/
-void  Mesh::ShowHistogram() const {
-	/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ShowHistogram)*/
-
-	const long kmax=10;
-	const double llmin = 0.5,llmax=2;
-	const double lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin);
-	long histo[kmax+1];
-	long i,it,k, nbedges =0;
-	for (i=0;i<=kmax;i++) histo[i]=0;
-	for (it=0;it<nbt;it++)
-	 if ( triangles[it].link) 
-		{
-
-		 for (int j=0;j<3;j++)
-			{
-			 Triangle *ta = triangles[it].TriangleAdj(j);
-			 if ( !ta || !ta->link || GetId(ta) >= it) 
-				{ 
-				 BamgVertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]];
-				 BamgVertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]];
-				 if ( !& vP || !&vQ) continue;
-				 R2 PQ = vQ.r - vP.r;
-				 double l = log(LengthInterpole(vP,vQ,PQ));
-				 nbedges++;
-				 k = (int) ((l - lmin)*delta);
-				 k = Min(Max(k,0L),kmax);
-				 histo[k]++;
-				}
-			}
-		}  
-	_printf_(" --- Histogram of the unit mesh,  nb of edges = " << nbedges << "\n");
-	_printf_("      length of edge in   | %% of edge  | Nb of edges \n"); 
-	_printf_("      --------------------+-------------+-------------\n"); 
-	for   (i=0;i<=kmax;i++){ 
-		if (i==0) _printf_( "      " << setw(10) << 0.);
-		else      _printf_( "      " << setw(10) << exp(lmin+i/delta));
-		if (i==kmax) _printf_("          +inf   ");
-		else      _printf_( "      " << setw(10) << exp(lmin+(i+1)/delta));
-		_printf_("|  " << setw(10) << (long((10000. * histo[i])/ nbedges)/100.) << " |\n");
-		_printf_("  " << histo[i] << "\n");
-	}
-	_printf_("      --------------------+-------------+-------------\n"); 
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/bamg/Mesh.h
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Mesh.h	(revision 16230)
+++ /issm/trunk-jpl/src/c/bamg/Mesh.h	(revision 16231)
@@ -83,7 +83,5 @@
 			long TriangleReferenceList(long*) const;
 			void TriangleIntNumbering(long* renumbering);
-			void ShowHistogram() const;
 			void CrackMesh(BamgOpts* bamgopts);
-			void ShowRegulaty() const;
 			void SmoothMetric(double raisonmax) ;
 			void BoundAnisotropy(double anisomax,double hminaniso= 1e-100) ;
@@ -96,5 +94,4 @@
 			long InsertNewPoints(long nbvold,long & NbTSwap) ; 
 			void TrianglesRenumberBySubDomain(bool justcompress=false);
-			void VerticesRenumber(long * renu);
 			void SmoothingVertex(int =3,double=0.3);
 			Metric MetricAt (const R2 &) const;
@@ -118,5 +115,4 @@
 			void BuildMetric1(BamgOpts* bamgopts);
 			void AddGeometryMetric(BamgOpts* bamgopts);
-			int  isCracked() const;
 			void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL);
 			void ReconstructExistingMesh();
Index: /issm/trunk-jpl/src/c/bamg/Triangle.cpp
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Triangle.cpp	(revision 16230)
+++ /issm/trunk-jpl/src/c/bamg/Triangle.cpp	(revision 16231)
@@ -331,9 +331,4 @@
 	void   Triangle::SetAllFlag(int a,int f){
 		AdjEdgeIndex[a] = (AdjEdgeIndex[a] &3) + (1020 & f);
-	}/*}}}*/
-	/*FUNCTION Triangle::SetDet{{{*/
-	void Triangle::SetDet() {
-		if(vertices[0] && vertices[1] && vertices[2])    det = bamg::det(*vertices[0],*vertices[1],*vertices[2]);
-		else det = -1; 
 	}/*}}}*/
 	/*FUNCTION Triangle::SetHidden{{{*/
Index: /issm/trunk-jpl/src/c/bamg/Triangle.h
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Triangle.h	(revision 16230)
+++ /issm/trunk-jpl/src/c/bamg/Triangle.h	(revision 16231)
@@ -63,5 +63,4 @@
 			void              SetMarkUnSwap(int a);
 			void              SetUnMarkUnSwap(int a);
-			void              SetDet();
 
 			//Inline methods
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16230)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16231)
@@ -3029,81 +3029,4 @@
 	/*clean-up*/
 	delete gauss;
-}
-/*}}}*/
-/*FUNCTION Penta::ReduceMatrixFS {{{*/
-void Penta::ReduceMatrixFS(IssmDouble* Ke_reduced, IssmDouble* Ke_temp){
-
-	int    i,j;
-	IssmDouble Kii[24][24];
-	IssmDouble Kib[24][3];
-	IssmDouble Kbb[3][3];
-	IssmDouble Kbi[3][24];
-	IssmDouble Kbbinv[3][3];
-	IssmDouble Kright[24][24];
-
-	/*Create the four matrices used for reduction */
-	for(i=0;i<24;i++){
-		for(j=0;j<24;j++){
-			Kii[i][j]=*(Ke_temp+27*i+j);
-		}
-		for(j=0;j<3;j++){
-			Kib[i][j]=*(Ke_temp+27*i+j+24);
-		}
-	}
-	for(i=0;i<3;i++){
-		for(j=0;j<24;j++){
-			Kbi[i][j]=*(Ke_temp+27*(i+24)+j);
-		}
-		for(j=0;j<3;j++){
-			Kbb[i][j]=*(Ke_temp+27*(i+24)+j+24);
-		}
-	}
-
-	/*Inverse the matrix corresponding to bubble part Kbb */
-	Matrix3x3Invert(&Kbbinv[0][0], &Kbb[0][0]);
-
-	/*Multiply matrices to create the reduce matrix Ke_reduced */
-	TripleMultiply(&Kib[0][0],24,3,0,
-				&Kbbinv[0][0],3,3,0,
-				&Kbi[0][0],3,24,0,
-				&Kright[0][0],0);
-
-	/*Affect value to the reduced matrix */
-	for(i=0;i<24;i++) for(j=0;j<24;j++) *(Ke_reduced+24*i+j)=Kii[i][j]-Kright[i][j];
-}
-/*}}}*/
-/*FUNCTION Penta::ReduceVectorFS {{{*/
-void Penta::ReduceVectorFS(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp){
-
-	int    i,j;
-	IssmDouble Pi[24];
-	IssmDouble Pb[3];
-	IssmDouble Kbb[3][3];
-	IssmDouble Kib[24][3];
-	IssmDouble Kbbinv[3][3];
-	IssmDouble Pright[24];
-
-	/*Create the four matrices used for reduction */
-	for(i=0;i<24;i++) Pi[i]=*(Pe_temp+i);
-	for(i=0;i<3;i++) Pb[i]=*(Pe_temp+i+24);
-	for(j=0;j<3;j++){
-		for(i=0;i<24;i++){
-			Kib[i][j]=*(Ke_temp+3*i+j);
-		}
-		for(i=0;i<3;i++){
-			Kbb[i][j]=*(Ke_temp+3*(i+24)+j);
-		}
-	}
-
-	/*Inverse the matrix corresponding to bubble part Kbb */
-	Matrix3x3Invert(&Kbbinv[0][0], &Kbb[0][0]);
-
-	/*Multiply matrices to create the reduce matrix Ke_reduced */
-	TripleMultiply(&Kib[0][0],24,3,0,
-				&Kbbinv[0][0],3,3,0,
-				&Pb[0],3,1,0,&Pright[0],0);
-
-	/*Affect value to the reduced matrix */
-	for(i=0;i<24;i++) *(Pe_reduced+i)=Pi[i]-Pright[i];
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16230)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16231)
@@ -224,6 +224,4 @@
 		bool           NoIceInElement(void); 
 		IssmDouble     MinEdgeLength(IssmDouble xyz_list[6][3]);
-		void	         ReduceMatrixFS(IssmDouble* Ke_reduced, IssmDouble* Ke_temp);
-		void	         ReduceVectorFS(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp);
 		void	         SetClone(int* minranks);
 		Tria*	         SpawnTria(int location);
Index: /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp	(revision 16230)
+++ /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp	(revision 16231)
@@ -733,21 +733,4 @@
 }
 /*}}}*/
-/*FUNCTION Riftfront::IsMaterialStable {{{*/
-int   Riftfront::IsMaterialStable(void){
-
-	IssmDouble converged=0;
-
-	this->inputs->GetInputValue(&converged,ConvergedEnum);
-
-	if(reCast<int,IssmDouble>(converged)){
-		/*ok, material non-linearity has converged. If that was already the case, we keep 
-		 * constraining the rift front. If it was not, and this is the first time the material 
-		 * has converged, we start constraining now!: */
-		this->material_converged=1;
-	}
-
-	return this->material_converged;
-}
-/*}}}*/
 /*FUNCTION Riftfront::MaxPenetration {{{*/
 int   Riftfront::MaxPenetration(IssmDouble* ppenetration){
Index: /issm/trunk-jpl/src/c/classes/Loads/Riftfront.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Riftfront.h	(revision 16230)
+++ /issm/trunk-jpl/src/c/classes/Loads/Riftfront.h	(revision 16231)
@@ -96,5 +96,4 @@
 		int   MaxPenetration(IssmDouble* ppenetration);
 		int   PotentialUnstableConstraint(int* punstable);
-		int   IsMaterialStable(void);
 		bool  IsInput(int name);
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h	(revision 16230)
+++ /issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h	(revision 16231)
@@ -21,11 +21,3 @@
 void RiftFreezeConstraints(Loads* loads,int analysis_type);
 
-/*rifts, trial and errors: */
-int  RiftIsPreStable(Loads* loads);
-void RiftSetPreStable(Loads* loads);
-void RiftPreConstrain(int* pnum_unstable_constraints,Loads* loads);
-void RiftMaxPenetrationInInputs(Loads* loads);
-int  RiftPotentialUnstableConstraints(Loads* loads);
-int  RiftIsMaterialStable(Loads* loads);
-
 #endif  /* _CONSTRAINTSSTATEX_H */
Index: /issm/trunk-jpl/src/c/modules/ConstraintsStatex/RiftConstraintsState.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ConstraintsStatex/RiftConstraintsState.cpp	(revision 16230)
+++ /issm/trunk-jpl/src/c/modules/ConstraintsStatex/RiftConstraintsState.cpp	(revision 16231)
@@ -159,189 +159,2 @@
 }
 /*}}}*/
-
-/*diverse trials and errors: */
-/*RiftIsMaterialStable(Loads* loads){{{*/
-int RiftIsMaterialStable(Loads* loads){
-
-	int i;
-
-	Riftfront* riftfront=NULL;
-	int found=0;
-	int mpi_found=0;
-
-	/*go though loads, and if non-linearity of the material has converged, let all penalties know: */
-	for (i=0;i<loads->Size();i++){
-
-		if(RiftfrontEnum==loads->GetEnum(i)){
-
-			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
-
-			if (riftfront->IsMaterialStable()){
-				found=1;
-				/*do not break! all penalties should get informed the non-linearity converged!*/
-			}
-		}
-	}
-
-	ISSM_MPI_Reduce (&found,&mpi_found,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&mpi_found,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
-	found=mpi_found;
-
-	return found;
-}
-/*}}}*/
-/*RiftIsPreStable(Loads* loads){{{*/
-int RiftIsPreStable(Loads* loads){
-
-	int i;
-
-	Riftfront* riftfront=NULL;
-	int found=0;
-	int mpi_found=0;
-
-	/*go though loads, and figure out if one of the penpair loads is still not stable: */
-	for (i=0;i<loads->Size();i++){
-
-		if(RiftfrontEnum==loads->GetEnum(i)){
-
-			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
-
-			if (riftfront->PreStable()==0){
-				found=1;
-				break;
-			}
-		}
-	}
-
-	ISSM_MPI_Reduce (&found,&mpi_found,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&mpi_found,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
-	found=mpi_found;
-
-	if (found){
-		/*We found an unstable constraint. : */
-		return 0;
-	}
-	else{
-		return 1;
-	}
-}
-/*}}}*/
-/*RiftSetPreStable(Loads* loads){{{*/
-void RiftSetPreStable(Loads* loads){
-
-	/*go though loads, and set loads to pre stable.:*/
-	for(int i=0;i<loads->Size();i++){
-		if(RiftfrontEnum==loads->GetEnum(i)){
-			Riftfront* riftfront=(Riftfront*)loads->GetObjectByOffset(i);
-			riftfront->SetPreStable();
-		}
-	}
-}
-/*}}}*/
-/*RiftPreConstrain(int* pnum_unstable_constraints,Loads* loads){{{*/
-void RiftPreConstrain(int* pnum_unstable_constraints,Loads* loads){
-
-	int			i;
-
-	/* generic object pointer: */
-	Riftfront* riftfront=NULL;
-
-	int unstable;
-	int sum_num_unstable_constraints;
-	int num_unstable_constraints=0;	
-
-	/*Enforce constraints: */
-	for (i=0;i<loads->Size();i++){
-
-		if (RiftfrontEnum==loads->GetEnum(i)){
-
-			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
-
-			riftfront->PreConstrain(&unstable);
-
-			num_unstable_constraints+=unstable;
-		}
-	}
-
-	ISSM_MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&sum_num_unstable_constraints,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
-	num_unstable_constraints=sum_num_unstable_constraints;
-
-	/*Assign output pointers: */
-	*pnum_unstable_constraints=num_unstable_constraints;
-
-}
-/*}}}*/
-/*RiftMaxPenetrationInInputs(Loads* loads){{{*/
-void RiftMaxPenetrationInInputs(Loads* loads){
-
-	int			i;
-
-	/* generic object pointer: */
-	Riftfront* riftfront=NULL;
-
-	/*rift penetration: */
-	IssmDouble max_penetration=0;
-	IssmDouble mpi_max_penetration;
-	IssmDouble penetration;
-
-	/*Ok, we are going to find the node pairs which are not penetrating, even though they 
-	 * are penalised. We will release only the one with has least <0 penetration. : */
-
-	max_penetration=0;
-	for (i=0;i<loads->Size();i++){
-
-		if (RiftfrontEnum==loads->GetEnum(i)){
-
-			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
-
-			riftfront->MaxPenetration(&penetration);
-
-			if (penetration>max_penetration)max_penetration=penetration;
-		}
-	}
-
-	ISSM_MPI_Reduce (&max_penetration,&mpi_max_penetration,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&mpi_max_penetration,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());                
-	max_penetration=mpi_max_penetration;
-
-	/*feed max_penetration to inputs: */
-	for(i=0;i<loads->Size();i++){
-		Load* load=(Load*)loads->GetObjectByOffset(i);
-		load->InputUpdateFromVector(&max_penetration,MaxPenetrationEnum,ConstantEnum);
-	}
-}
-/*}}}*/
-/*RiftPotentialUnstableConstraints(Loads* loads){{{*/
-int RiftPotentialUnstableConstraints(Loads* loads){
-
-	int			i;
-
-	/* generic object pointer: */
-	Riftfront* riftfront=NULL;
-
-	/*Ok, we are going to find the node pairs which are not penetrating, even though they 
-	 * are penalised. We will release only the one with has least <0 penetration. : */
-	int unstable=0;
-	int sum_num_unstable_constraints=0;
-	int num_unstable_constraints=0;
-
-	for (i=0;i<loads->Size();i++){
-
-		if (RiftfrontEnum==loads->GetEnum(i)){
-
-			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
-
-			riftfront->PotentialUnstableConstraint(&unstable);
-
-			num_unstable_constraints+=unstable;
-		}
-	}
-
-	ISSM_MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&sum_num_unstable_constraints,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
-	num_unstable_constraints=sum_num_unstable_constraints;
-
-	return num_unstable_constraints;
-}
-/*}}}*/
Index: /issm/trunk-jpl/src/c/shared/TriMesh/TriMeshUtils.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/TriMesh/TriMeshUtils.cpp	(revision 16230)
+++ /issm/trunk-jpl/src/c/shared/TriMesh/TriMeshUtils.cpp	(revision 16231)
@@ -490,178 +490,4 @@
 	return noerr;
 }/*}}}*/
-/*FUNCTION RemoveRifts{{{*/
-int RemoveRifts(int** pindex,double** px,double** py,int* pnods,int** psegments,int* pnumsegs,int numrifts1,int* rifts1numsegs,int** rifts1segments,int** rifts1pairs,int nel){
-
-	int noerr=1;
-	int i,j,k,counter,counter1,counter2;
-
-	/*intermediary: */
-	int    *riftsegments = NULL;
-	int    *riftpairs    = NULL;
-	int     node1,node2,node3,node4,temp_node;
-	int     el2;
-	int     newnods; //temporary # node counter.
-	double  xmin,ymin;
-	double *xreal        = NULL;
-	double *yreal        = NULL;
-	int    *nodes        = NULL;
-	int    *mergingnodes = NULL;
-	int     max_size;
-	int     redundant;
-
-	/*Recover input: */
-	int    *index    = *pindex;
-	double *x        = *px;
-	double *y        = *py;
-	int     nods     = *pnods;     ;
-	int    *segments = *psegments;
-	int     numsegs  = *pnumsegs;
-
-	/*initialize newnods : */
-	newnods=nods;
-
-	/*Figure out a unique value to flag x and y for node removal: */
-	xmin=x[0];
-	ymin=y[0];
-	for (i=0;i<nods;i++){
-		if (x[i]<xmin)xmin=x[i];
-		if (y[i]<ymin)ymin=y[i];
-	}
-	xmin=xmin-fabs(xmin); 
-	ymin=ymin-fabs(ymin);
-
-	/*Initialize two arrays, one for nodes that are going to be merged, the other with corresponding nodes being merge into: */
-	max_size=0;
-	for (i=0;i<numrifts1;i++){
-		max_size+=rifts1numsegs[i];
-	}
-	nodes=xNew<int>(max_size);
-	mergingnodes=xNew<int>(max_size);
-
-	/*Go through the rifts segments, and identify which node we are going to merge with its counterpart on the other side 
-	 *of the rift. The way we identify this node is by looking at the element pairs, and the corresponding nodes: */
-	counter=0;
-	for (i=0;i<numrifts1;i++){
-		riftsegments=rifts1segments[i];
-		riftpairs=rifts1pairs[i];
-		for (j=0;j<rifts1numsegs[i];j++){
-			el2=riftpairs[2*j+1];
-			node1=(int)*(riftsegments+3*j+0);
-			node2=(int)*(riftsegments+3*j+1);
-			/*Summary, el1 and el2 are facing one another across the rift. node1 and node2 belong to el1 and 
-			 *are located on the rift. Find node3 and node4, nodes belonging to el2 and located on the rift: */
-			for (k=0;k<rifts1numsegs[i];k++){
-				if (*(riftsegments+3*k+2)==el2){
-					node3=*(riftsegments+3*k+0);
-					node4=*(riftsegments+3*k+1);
-					break;
-				}
-			}
-			/* Make sure node3 faces node1 and node4 faces node2: */
-			if ((x[node1]==x[node4]) && (y[node1]==y[node4])){
-				/*Swap node3 and node4:*/
-				temp_node=node3;
-				node3=node4;
-				node4=temp_node;
-			}
-			/* Is any of these two node pairs on the tip of a rift, in which case, we don't include it in nodes or mergingnodes: */
-			if ((node1==node3) || (node2==node4)){
-				if(node1!=node3){
-					/*Add node1 and node3 to nodes and mergingnodes if they have not already been added: */
-					redundant=0;
-					for (k=0;k<counter;k++){
-						if ((mergingnodes[k]==node1) || (nodes[k]==node1))redundant=1;
-					}
-					if(!redundant){
-						/*Ok, add node1 to nodes, and node3 to mergingnodes: */
-						nodes[counter]=node1;
-						mergingnodes[counter]=node3;
-						counter++;
-					}
-				}
-				if(node2!=node4){
-					/*Add node2 and node4 to nodes and mergingnodes if they have not already been added: */
-					redundant=0;
-					for (k=0;k<counter;k++){
-						if ((mergingnodes[k]==node2) || (nodes[k]==node2))redundant=1;
-					}
-					if(!redundant){
-						/*Ok, add node2 to nodes, and node4 to mergingnodes: */
-						nodes[counter]=node2;
-						mergingnodes[counter]=node4;
-						counter++;
-					}
-				}
-			}
-			else{
-				/*Check that node1 is not already present in the mergingnodes: */
-				redundant=0;
-				for (k=0;k<counter;k++){
-					if ((mergingnodes[k]==node1) || (nodes[k]==node1))redundant=1;
-				}
-				if(!redundant){
-					/*Ok, add node1 to nodes, and node3 to mergingnodes: */
-					nodes[counter]=node1;
-					mergingnodes[counter]=node3;
-					counter++;
-				}
-				/*Check that node2 is not already present in the mergingnodes: */
-				redundant=0;
-				for (k=0;k<counter;k++){
-					if ((mergingnodes[k]==node1) || (nodes[k]==node1))redundant=1;
-				}
-				if(!redundant){
-					/*Ok, add node2 to nodes, and node4 to mergingnodes: */
-					nodes[counter]=node2;
-					mergingnodes[counter]=node4;
-					counter++;
-				}
-			}
-		}
-	}
-
-	/*Ok, we have counter pairs of nodes (nodes and mergingnodes): start merging nodes in the triangulation: */
-	newnods=nods;
-	for (i=0;i<counter;i++){
-		node1=nodes[i];
-		node3=mergingnodes[i];
-		/*Merge node3 into node1: */ 
-		x[node3]=xmin; //flag  for later removal from x
-		y[node3]=ymin; //flag  for later removal from y
-		newnods--;
-		for (k=0;k<nel;k++){
-			if (*(index+3*k+0)==node3)*(index+3*k+0)=node1;
-			if (*(index+3*k+1)==node3)*(index+3*k+1)=node1;
-			if (*(index+3*k+2)==node3)*(index+3*k+2)=node1;
-		}
-	}
-
-	/*Reallocate x and y: */
-	xreal=xReNew<double>(x,nods,newnods);
-	yreal=xReNew<double>(y,nods,newnods);
-	counter1=0;
-	counter2=0;
-	for (i=0;i<nods;i++){
-		if (x[i]!=xmin){
-			xreal[counter1]=x[i];
-			counter1++;
-		}
-		if (y[i]!=ymin){
-			yreal[counter2]=y[i];
-			counter2++;
-		}
-	}
-	xDelete<double>(x); x=xreal;
-	xDelete<double>(y); y=yreal;
-
-	/*Assign output pointers:*/
-	*pindex=index;
-	*px=x;
-	*py=y;
-	*pnods=newnods;
-	*psegments=segments;
-	*pnumsegs=numsegs;
-	return noerr;
-}/*}}}*/
 /*FUNCTION IsRiftPresent{{{*/
 int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs){
Index: /issm/trunk-jpl/src/c/shared/TriMesh/trimesh.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/TriMesh/trimesh.h	(revision 16230)
+++ /issm/trunk-jpl/src/c/shared/TriMesh/trimesh.h	(revision 16231)
@@ -23,5 +23,4 @@
 int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel);
 int FindElement(double A,double B,int* index,int nel);
-int RemoveRifts(int** pindex,double** px,double** py,int* pnods,int** psegments,int* pnumsegs,int numrifts1,int* rifts1numsegs,int** rifts1segments,double** rifts1pairs,int nel);
 int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs);
 int SplitRiftSegments(int** psegments,int** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,int*** priftssegments,int numrifts,int nods,int nels);
Index: sm/trunk-jpl/src/c/toolkits/petsc/patches/PetscMatrixToDoubleMatrix.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/patches/PetscMatrixToDoubleMatrix.cpp	(revision 16230)
+++ 	(revision )
@@ -1,47 +1,0 @@
-/* \file PetscMatrixToDoubleMatrix.cpp
- * \brief: convert a sparse or dense Petsc matrix into a matlab matrix
- */
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#include <string>
-
-/*Petsc includes: */
-#include <petscmat.h>
-#include <petscvec.h>
-#include <petscksp.h>
-
-/*Petsc includes: */
-#include "../../../shared/shared.h"
-
-void PetscMatrixToDoubleMatrix(double** pmatrix, int* prows, int* pcols,Mat petsc_matrix){
-
-	/*output: */
-	int     i;
-	double* matrix=NULL;
-	int     rows,cols;
-
-	/*Some needed information: */
-	MatGetSize(petsc_matrix,&rows,&cols);
-
-	int* idxm=xNew<int>(rows);
-	int* idxn=xNew<int>(cols);
-
-	for(i=0;i<rows;i++)idxm[i]=i;
-	for(i=0;i<cols;i++)idxn[i]=i;
-
-	matrix=xNew<double>(rows*cols);
-	MatGetValues(petsc_matrix,rows,idxm,cols,idxn,matrix);
-
-	xDelete<int>(idxm);
-	xDelete<int>(idxn);
-
-	/*Assign output pointers: */
-	*pmatrix=matrix;
-	*prows=rows;
-	*pcols=cols;
-}
Index: sm/trunk-jpl/src/c/toolkits/petsc/patches/PetscVectorToDoubleVector.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/patches/PetscVectorToDoubleVector.cpp	(revision 16230)
+++ 	(revision )
@@ -1,42 +1,0 @@
-/* \file PetscVectorToDoubleVector.cpp
- */
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#include <string>
-
-/*Petsc includes: */
-#include <petscmat.h>
-#include <petscvec.h>
-#include <petscksp.h>
-
-#include "../../../shared/shared.h"
-
-void PetscVectorToDoubleVector(double** pvector, int* prows, Vec petsc_vector){
-
-	int     rows;
-	double *vector = NULL;
-
-	/*Get size of vector: */
-	if(petsc_vector){
-		VecGetSize(petsc_vector,&rows);
-		if(rows){
-			int* idxm=xNew<int>(rows);
-			vector=xNew<double>(rows);
-			for(int i=0;i<rows;i++)idxm[i]=i;
-			VecGetValues(petsc_vector,rows,idxm,vector);
-			xDelete<int>(idxm);
-		}
-	}
-	else{
-		rows=0;
-	}
-
-	/*Assign output pointers: */
-	*pvector=vector;
-	*prows=rows;
-}
Index: /issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h	(revision 16230)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h	(revision 16231)
@@ -39,6 +39,3 @@
 MatType ISSMToPetscMatrixType(MatrixType type);
 
-void PetscMatrixToDoubleMatrix(double** pmatrix, int* prows, int* pcols,Mat matrix);
-void PetscVectorToDoubleVector(double** pvector, int* prows, Vec vector);
-
 #endif