Index: /issm/trunk-jpl/src/c/CMakeLists.txt
===================================================================
--- /issm/trunk-jpl/src/c/CMakeLists.txt	(revision 16227)
+++ /issm/trunk-jpl/src/c/CMakeLists.txt	(revision 16228)
@@ -517,5 +517,4 @@
 					./toolkits/petsc/patches/PetscMatrixToDoubleMatrix.cpp
 					./toolkits/petsc/patches/PetscVectorToDoubleVector.cpp
-					./toolkits/petsc/patches/VecDuplicatePatch.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 16227)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 16228)
@@ -197,5 +197,4 @@
 					./shared/Exceptions/exceptions.h\
 					./shared/Exceptions/Exceptions.cpp\
-					./shared/Exceptions/exprintf.cpp\
 					./shared/Sorting/binary_search.cpp\
 					./shared/Sorting/sorting.h\
@@ -770,5 +769,4 @@
 					./toolkits/petsc/patches/PetscMatrixToDoubleMatrix.cpp\
 					./toolkits/petsc/patches/PetscVectorToDoubleVector.cpp\
-					./toolkits/petsc/patches/VecDuplicatePatch.cpp\
 					./toolkits/petsc/patches/KSPFree.cpp\
 					./toolkits/petsc/patches/MatFree.cpp\
Index: /issm/trunk-jpl/src/c/bamg/Mesh.cpp
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Mesh.cpp	(revision 16227)
+++ /issm/trunk-jpl/src/c/bamg/Mesh.cpp	(revision 16228)
@@ -4182,510 +4182,4 @@
 }
 /*}}}*/
-	/*FUNCTION Mesh::SplitElement{{{*/
-	int  Mesh::SplitElement(int choice){
-		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/SplitElement)*/
-
-		long int verbose=0;
-
-		Direction NoDirOfSearch;
-		const  int withBackground = &BTh != this && &BTh;
-
-		TrianglesRenumberBySubDomain();
-		//int nswap =0;
-		const long nfortria( choice ? 4 : 6);
-		if(withBackground) 
-		  {
-			BTh.SetVertexFieldOn();
-			SetVertexFieldOnBTh();
-		  }
-		else
-		 BTh.SetVertexFieldOn();
-
-		long newnbt=0,newnbv=0;
-		long * kedge = 0;
-		long newnbq=nbq;
-		long * ksplit= 0, * ksplitarray=0;
-		long kkk=0;
-		int ret =0;
-		if (maxnbv<nbv+nbe) return 1;//   
-		// 1) create  the new points by spliting the internal edges 
-		// set the 
-		long nbvold = nbv;
-		long nbtold = nbt;
-		long nbtoutold  = nbtout;
-		long  NbEdgeOnGeom=0;
-		long i;
-
-		nbt = nbt - nbtout; // remove all the  the ouside triangles 
-		long nbtsave = nbt;
-		Triangle * lastT = triangles + nbt;
-		for (i=0;i<nbe;i++)
-		 if(edges[i].GeomEdgeHook) NbEdgeOnGeom++;
-		long newnbe=nbe+nbe;
-		//  long newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex;
-		long newNbVerticesOnGeomEdge=NbVerticesOnGeomEdge+NbEdgeOnGeom;
-		// long newNbVertexOnBThVertex=NbVertexOnBThVertex;
-		long newNbVertexOnBThEdge=withBackground ? NbVertexOnBThEdge+NbEdgeOnGeom:0;
-
-		// do allocation for pointeur to the geometry and background
-		VertexOnGeom * newVerticesOnGeomEdge = new VertexOnGeom[newNbVerticesOnGeomEdge];
-		VertexOnEdge *newVertexOnBThEdge = newNbVertexOnBThEdge ?  new VertexOnEdge[newNbVertexOnBThEdge]:0;
-		if (NbVerticesOnGeomEdge)
-		 memcpy(newVerticesOnGeomEdge,VerticesOnGeomEdge,sizeof(VertexOnGeom)*NbVerticesOnGeomEdge);
-		if (NbVertexOnBThEdge)
-		 memcpy(newVertexOnBThEdge,VertexOnBThEdge,sizeof(VertexOnEdge)*NbVertexOnBThEdge);
-		Edge *newedges = new Edge [newnbe];
-		//  memcpy(newedges,edges,sizeof(Edge)*nbe);
-		SetOfEdges4 * edge4= new SetOfEdges4(nbe,nbv);
-		long k=nbv;
-		long kk=0;
-		long kvb = NbVertexOnBThEdge;
-		long kvg = NbVerticesOnGeomEdge;
-		long ie =0;
-		Edge ** edgesGtoB=0;
-		if (withBackground)
-		 edgesGtoB= BTh.MakeGeomEdgeToEdge();
-		long ferr=0;
-		for (i=0;i<nbe;i++)
-		 newedges[ie].GeomEdgeHook=0;
-
-		for (i=0;i<nbe;i++)
-		  {
-			GeomEdge *ong =  edges[i].GeomEdgeHook;
-
-			newedges[ie]=edges[i];
-			newedges[ie].adj[0]=newedges+(edges[i].adj[0]-edges) ;
-			newedges[ie].adj[1]=newedges + ie +1;
-			R2 A = edges[i][0],B = edges[i][1];
-
-			kk += (i == edge4->SortAndAdd(GetId(edges[i][0]),GetId(edges[i][1])));
-			if (ong) // a geometrical edges 
-			  { 
-				if (withBackground){
-					// walk on back ground mesh 
-					//  newVertexOnBThEdge[ibe++] = VertexOnEdge(vertices[k],bedge,absicsseonBedge); 
-					// a faire -- difficile 
-					// the first PB is to now a background edge between the 2 vertices
-					if (!edgesGtoB){
-						_error_("!edgesGtoB");
-					}
-					ong= ProjectOnCurve(*edgesGtoB[Gh.GetId(edges[i].GeomEdgeHook)],
-								edges[i][0],edges[i][1],0.5,vertices[k],
-								newVertexOnBThEdge[kvb],
-								newVerticesOnGeomEdge[kvg++]);
-					vertices[k].ReferenceNumber= edges[i].ReferenceNumber;
-					vertices[k].DirOfSearch =   NoDirOfSearch;        
-					;
-					// get the Info on background mesh 
-					double s =        newVertexOnBThEdge[kvb];
-					BamgVertex &  bv0  = newVertexOnBThEdge[kvb][0];
-					BamgVertex &  bv1  = newVertexOnBThEdge[kvb][1];
-					// compute the metrix of the new points 
-					vertices[k].m =  Metric(1-s,bv0,s,bv1); 
-					kvb++;
-				  }
-				else 
-				  {
-					ong=Gh.ProjectOnCurve(edges[i],
-								0.5,vertices[k],newVerticesOnGeomEdge[kvg++]);
-					// vertices[k].i = R2ToI2( vertices[k].r);
-					vertices[k].ReferenceNumber = edges[i].ReferenceNumber;
-					vertices[k].DirOfSearch = NoDirOfSearch;
-					vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);	      
-				  }  
-			  }
-			else // straigth line edge ---
-			  { 
-				vertices[k].r = ((R2) edges[i][0] + (R2)  edges[i][1] )*0.5;
-				vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);
-				vertices[k].GeomEdgeHook = 0;
-			  }
-			//vertices[k].i = R2ToI2( vertices[k].r);
-			R2 AB =  vertices[k].r;
-			R2 AA = (A+AB)*0.5;
-			R2 BB = (AB+B)*0.5;
-			vertices[k].ReferenceNumber = edges[i].ReferenceNumber;
-			vertices[k].DirOfSearch = NoDirOfSearch;
-
-			newedges[ie].GeomEdgeHook = Gh.Containing(AA,ong);
-			newedges[ie++].v[1]=vertices+k;
-
-			newedges[ie]=edges[i];
-			newedges[ie].adj[0]=newedges + ie -1;
-			newedges[ie].adj[1]=newedges+(edges[i].adj[1]-edges) ;
-			newedges[ie].GeomEdgeHook =  Gh.Containing(BB,ong);
-			newedges[ie++].v[0]=vertices+k;
-			k++;
-		  }
-		if (edgesGtoB) delete [] edgesGtoB;
-		edgesGtoB=0;
-
-		newnbv=k;
-		newNbVerticesOnGeomEdge=kvg;
-		if (newnbv> maxnbv) goto Error;// bug 
-
-		nbv = k;
-
-		kedge = new long[3*nbt+1];
-		ksplitarray = new long[nbt+1];
-		ksplit = ksplitarray +1; // because ksplit[-1] == ksplitarray[0]
-
-		for (i=0;i<3*nbt;i++)
-		 kedge[i]=-1;
-
-		//  
-
-		for (i=0;i<nbt;i++) {
-			Triangle & t = triangles[i];
-			if (!t.link){
-				_error_("!t.link");
-			}
-			for(int j=0;j<3;j++)
-			  {
-				const AdjacentTriangle ta = t.Adj(j);
-				const Triangle & tt = ta;
-				if (&tt >= lastT)
-				 t.SetAdj2(j,0,0);// unset adj
-				const BamgVertex & v0 = t[VerticesOfTriangularEdge[j][0]];
-				const BamgVertex & v1 = t[VerticesOfTriangularEdge[j][1]];
-				long  ke =edge4->SortAndFind(GetId(v0),GetId(v1));
-				if (ke>0) 
-				  {
-					long ii = GetId(tt);
-					int  jj = ta;
-					long ks = ke + nbvold;
-					kedge[3*i+j] = ks;
-					if (ii<nbt) // good triangle
-					 kedge[3*ii+jj] = ks;
-					BamgVertex &A=vertices[ks];
-					double aa,bb,cc,dd;
-					if ((dd=Area2(v0.r,v1.r,A.r)) >=0){
-						// warning PB roundoff error 
-						if (t.link && ( (aa=Area2( A.r    , t[1].r , t[2].r )) < 0.0 
-										||   (bb=Area2( t[0].r , A.r    , t[2].r )) < 0.0  
-										||   (cc=Area2( t[0].r , t[1].r , A.r    )) < 0.0)){
-							_printf_(ke + nbvold << " not in triangle " << i << " In= " << !!t.link << " " << aa << " " << bb << " " << cc << " " << dd << "\n");
-							_error_("Number of triangles with P2 interpolation Problem");
-						}
-					}
-					else {
-						if (tt.link && ( (aa=Area2( A.r     , tt[1].r , tt[2].r )) < 0 
-										||   (bb=Area2( tt[0].r , A.r     , tt[2].r )) < 0 
-										||   (cc=Area2( tt[0].r , tt[1].r , A.r     )) < 0)){
-							_printf_(ke + nbvold << " not in triangle " << ii << " In= " << !!tt.link << " " << aa << " " << bb << " " << cc << " " << dd << "\n");
-							_error_("Number of triangles with P2 interpolation Problem");
-						}
-					} 
-				  }
-			  }
-		}
-
-		for (i=0;i<nbt;i++){
-			ksplit[i]=1; // no split by default
-			const Triangle & t = triangles[ i];
-			int nbsplitedge =0;
-			int nbinvisible =0;
-			int invisibleedge=0;
-			int kkk[3];      
-			for (int j=0;j<3;j++)
-			  {
-				if (t.Hidden(j)) invisibleedge=j,nbinvisible++;
-
-				const AdjacentTriangle ta = t.Adj(j);
-				const Triangle & tt = ta;
-
-				const BamgVertex & v0 = t[VerticesOfTriangularEdge[j][0]];
-				const BamgVertex & v1 = t[VerticesOfTriangularEdge[j][1]];
-				if ( kedge[3*i+j] < 0) 
-				  {
-					long  ke =edge4->SortAndFind(GetId(v0),GetId(v1));
-					if (ke<0) // new 
-					  {
-						if (&tt) // internal triangles all the boundary 
-						  { // new internal edges 
-							long ii = GetId(tt);
-							int  jj = ta;
-
-							kedge[3*i+j]=k;// save the vertex number 
-							kedge[3*ii+jj]=k;
-							if (k<maxnbv) 
-							  {
-								vertices[k].r = ((R2) v0+(R2) v1 )/2;
-								//vertices[k].i = R2ToI2( vertices[k].r);
-								vertices[k].ReferenceNumber=0;
-								vertices[k].DirOfSearch =NoDirOfSearch;
-								vertices[k].m =  Metric(0.5,v0,0.5,v1);
-							  }
-							k++;
-							kkk[nbsplitedge++]=j;		      
-						  } // tt 
-						else
-						 _error_("Bug...");
-					  } // ke<0	       
-					else
-					  { // ke >=0
-						kedge[3*i+j]=nbvold+ke;
-						kkk[nbsplitedge++]=j;// previously splited
-					  }
-				  }
-				else 
-				 kkk[nbsplitedge++]=j;// previously splited
-
-			  } 
-			if (nbinvisible>=2){
-				_error_("nbinvisible>=2");
-			}
-			switch (nbsplitedge) {
-				case 0: ksplit[i]=10; newnbt++; break;   // nosplit
-				case 1: ksplit[i]=20+kkk[0];newnbt += 2; break; // split in 2 
-				case 2: ksplit[i]=30+3-kkk[0]-kkk[1];newnbt += 3; break; // split in 3 
-				case 3:
-						  if (nbinvisible) ksplit[i]=40+invisibleedge,newnbt += 4;
-						  else   ksplit[i]=10*nfortria,newnbt+=nfortria;
-						  break;
-			} 
-			if (ksplit[i]<40){
-				_error_("ksplit[i]<40");
-			}
-		  }
-		//  now do the element split
-		newnbq = 4*nbq;
-		nbv = k;
-		kkk = nbt;
-		ksplit[-1] = nbt;
-		// look on  old true  triangles 
-
-		for (i=0;i<nbtsave;i++){
-			int  nbmkadj=0;
-			long mkadj [100];
-			mkadj[0]=i;
-			long kk=ksplit[i]/10;
-			int  ke=(int) (ksplit[i]%10);
-			if (kk>=7 || kk<=0){
-				_error_("kk>=7 || kk<=0");
-			}
-
-			// def the numbering   k (edge) i vertex 
-			int k0 = ke;
-			int k1 = NextEdge[k0];
-			int k2 = PreviousEdge[k0];
-			int i0 = OppositeVertex[k0];
-			int i1 = OppositeVertex[k1];
-			int i2 = OppositeVertex[k2];
-
-			Triangle &t0=triangles[i];
-			BamgVertex * v0=t0(i0);           
-			BamgVertex * v1=t0(i1);           
-			BamgVertex * v2=t0(i2);
-
-			if (nbmkadj>=10){
-				_error_("nbmkadj>=10");
-			}
-			// --------------------------
-			AdjacentTriangle ta0(t0.Adj(i0)),ta1(t0.Adj(i1)),ta2(t0.Adj(i2));
-			// save the flag Hidden
-			int hid[]={t0.Hidden(0),t0.Hidden(1),t0.Hidden(2)};
-			// un set all adj -- save Hidden flag --
-			t0.SetAdj2(0,0,hid[0]);
-			t0.SetAdj2(1,0,hid[1]);
-			t0.SetAdj2(2,0,hid[2]);
-			// --  remake 
-			switch  (kk) {
-				case 1: break;// nothing 
-				case 2: // 
-						  {
-							Triangle &t1=triangles[kkk++];
-							t1=t0;
-							if (kedge[3*i+i0]<0){
-								_error_("kedge[3*i+i0]<0");
-							}
-							BamgVertex * v3 = vertices + kedge[3*i+k0];
-
-							t0(i2) = v3;
-							t1(i1) = v3;
-							t0.SetAllFlag(k2,0);
-							t1.SetAllFlag(k1,0);
-						  } 
-						break; 
-				case 3: //
-						  {
-							Triangle &t1=triangles[kkk++];
-							Triangle &t2=triangles[kkk++];
-							t2=t1=t0;
-							if (kedge[3*i+k1]<0){
-								_error_("kedge[3*i+k1]<0");
-							}
-							if (kedge[3*i+k2]<0){
-								_error_("kedge[3*i+k2]<0");
-							}
-
-							BamgVertex * v01 = vertices + kedge[3*i+k2];
-							BamgVertex * v02 = vertices + kedge[3*i+k1]; 
-							t0(i1) = v01; 
-							t0(i2) = v02; 
-							t1(i2) = v02;
-							t1(i0) = v01; 
-							t2(i0) = v02; 
-							t0.SetAllFlag(k0,0);
-							t1.SetAllFlag(k1,0);
-							t1.SetAllFlag(k0,0);
-							t2.SetAllFlag(k2,0);
-						  } 
-						break;
-				case 4: // 
-				case 6: // split in 4 
-						  {
-							Triangle &t1=triangles[kkk++];
-							Triangle &t2=triangles[kkk++];
-							Triangle &t3=triangles[kkk++];
-							t3=t2=t1=t0;
-							if (kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0){
-								_error_("kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0");
-							}
-							BamgVertex * v12 = vertices + kedge[3*i+k0];
-							BamgVertex * v02 = vertices + kedge[3*i+k1]; 
-							BamgVertex * v01 = vertices + kedge[3*i+k2];
-							t0(i1) = v01;
-							t0(i2) = v02;
-							t0.SetAllFlag(k0,hid[k0]);
-
-							t1(i0) = v01;
-							t1(i2) = v12;
-							t0.SetAllFlag(k1,hid[k1]);
-
-							t2(i0) = v02;
-							t2(i1) = v12;
-							t2.SetAllFlag(k2,hid[k2]);
-
-							t3(i0) = v12;
-							t3(i1) = v02;
-							t3(i2) = v01;
-
-							t3.SetAllFlag(0,hid[0]);	   
-							t3.SetAllFlag(1,hid[1]);	   
-							t3.SetAllFlag(2,hid[2]);
-
-							if ( kk == 6)
-							  {
-
-								Triangle &t4=triangles[kkk++];
-								Triangle &t5=triangles[kkk++];
-
-								t4 = t3;
-								t5 = t3;
-
-								t0.SetHidden(k0);
-								t1.SetHidden(k1);
-								t2.SetHidden(k2);
-								t3.SetHidden(0);
-								t4.SetHidden(1);
-								t5.SetHidden(2);
-
-								if (nbv < maxnbv ) 
-								  {
-									vertices[nbv].r = ((R2) *v01 + (R2) *v12  + (R2) *v02 ) / 3.0;
-									vertices[nbv].ReferenceNumber =0;
-									vertices[nbv].DirOfSearch =NoDirOfSearch;
-									//vertices[nbv].i = R2ToI2(vertices[nbv].r);
-									double a3[]={1./3.,1./3.,1./3.};
-									vertices[nbv].m = Metric(a3,v0->m,v1->m,v2->m);
-									BamgVertex * vc =  vertices +nbv++;
-									t3(i0) = vc;
-									t4(i1) = vc;
-									t5(i2) = vc;
-
-								  }
-								else
-								 goto Error; 
-							  }
-
-						  } 
-						break;         
-			}
-
-			// save all the new triangles
-			mkadj[nbmkadj++]=i;
-			long jj;
-			if (t0.link) 
-			 for (jj=nbt;jj<kkk;jj++)
-				{
-				 triangles[jj].link=t0.link;
-				 t0.link= triangles+jj;
-				 mkadj[nbmkadj++]=jj;
-				}
-			if (nbmkadj>13){// 13 = 6 + 4 +
-				_error_("nbmkadj>13");
-			}
-
-			if (kk==6)  newnbq+=3;
-			for (jj=ksplit[i-1];jj<kkk;jj++) nbt = kkk;
-			ksplit[i]= nbt; // save last adresse of the new triangles
-			kkk = nbt;
-		  }
-
-		for (i=0;i<nbv;i++) vertices[i].m = vertices[i].m*2.;
-
-		if(withBackground)
-		 for (i=0;i<BTh.nbv;i++)
-		  BTh.vertices[i].m =  BTh.vertices[i].m*2.;
-
-		ret = 2;
-		if (nbt>= maxnbt) goto Error; // bug 
-		if (nbv>= maxnbv) goto Error; // bug 
-		// generation of the new triangles 
-
-		SetIntCoor("In SplitElement"); 
-
-		CreateSingleVertexToTriangleConnectivity();
-		if(withBackground)
-		 BTh.CreateSingleVertexToTriangleConnectivity();
-
-		delete [] edges;
-		edges = newedges;
-		nbe = newnbe;
-		nbq = newnbq;
-
-		for (i=0;i<nbsubdomains;i++)
-		  { 
-			long k = subdomains[i].edge- edges;
-			subdomains[i].edge =  edges+2*k; // spilt all edge in 2 
-		  }
-
-		if (ksplitarray) delete [] ksplitarray;
-		if (kedge) delete [] kedge;
-		if (edge4) delete edge4;
-		if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge;
-		VerticesOnGeomEdge= newVerticesOnGeomEdge;
-		if(VertexOnBThEdge) delete []  VertexOnBThEdge;
-		VertexOnBThEdge = newVertexOnBThEdge;
-		NbVerticesOnGeomEdge = newNbVerticesOnGeomEdge;
-		NbVertexOnBThEdge=newNbVertexOnBThEdge;
-		//  CreateSingleVertexToTriangleConnectivity();
-
-		ReconstructExistingMesh();
-
-		if (verbose>2){
-			_printf_("   number of quadrilaterals    = " << nbq << "\n");
-			_printf_("   number of triangles         = " << nbt-nbtout- nbq*2 << "\n");
-			_printf_("   number of outside triangles = " << nbtout << "\n");
-		}
-
-		return 0; //ok
-
-Error:
-		nbv = nbvold;
-		nbt = nbtold;
-		nbtout = nbtoutold;
-		// cleaning memory ---
-		delete [] newedges;
-		if (ksplitarray) delete [] ksplitarray;
-		if (kedge) delete [] kedge;
-		if (newVerticesOnGeomEdge) delete [] newVerticesOnGeomEdge;
-		if (edge4) delete edge4;
-		if(newVertexOnBThEdge) delete []  newVertexOnBThEdge;
-
-		return ret; // ok 
-	}
-	/*}}}*/
 /*FUNCTION Mesh::SplitInternalEdgeWithBorderVertices{{{*/
 long  Mesh::SplitInternalEdgeWithBorderVertices(){
Index: /issm/trunk-jpl/src/c/bamg/Mesh.h
===================================================================
--- /issm/trunk-jpl/src/c/bamg/Mesh.h	(revision 16227)
+++ /issm/trunk-jpl/src/c/bamg/Mesh.h	(revision 16228)
@@ -92,5 +92,4 @@
 			long SplitInternalEdgeWithBorderVertices();
 			void MakeQuadrangles(double costheta);
-			int  SplitElement(int choice);
 			void MakeBamgQuadtree();
 			void NewPoints(Mesh &,BamgOpts* bamgopts,int KeepVertices=1);
Index: /issm/trunk-jpl/src/c/bamg/SetOfE4.cpp
===================================================================
--- /issm/trunk-jpl/src/c/bamg/SetOfE4.cpp	(revision 16227)
+++ /issm/trunk-jpl/src/c/bamg/SetOfE4.cpp	(revision 16228)
@@ -103,9 +103,4 @@
 	}
 	/*}}}*/
-	/*FUNCTION  SetOfEdges4::newarete{{{*/
-	long SetOfEdges4::newarete(long k){
-		return NbOfEdges == k+1;
-	}
-	/*}}}*/
 	/*FUNCTION  SetOfEdges4::SortAndAdd{{{*/
 	long SetOfEdges4::SortAndAdd (long ii,long jj) {
Index: /issm/trunk-jpl/src/c/bamg/SetOfE4.h
===================================================================
--- /issm/trunk-jpl/src/c/bamg/SetOfE4.h	(revision 16227)
+++ /issm/trunk-jpl/src/c/bamg/SetOfE4.h	(revision 16228)
@@ -37,5 +37,4 @@
 			long i(long k);
 			long j(long k);
-			long newarete(long k);
 	};
 }
Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 16227)
+++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp	(revision 16228)
@@ -820,7 +820,2 @@
 }
 /*}}}*/
-/*FUNCTION Pengrid::UpdateInputs {{{*/
-void  Pengrid::UpdateInputs(IssmDouble* solution){
-	_error_("not supported yet!");
-}
-/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h	(revision 16227)
+++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h	(revision 16228)
@@ -103,5 +103,4 @@
 		#endif
 		void  ConstraintActivate(int* punstable);
-		void  UpdateInputs(IssmDouble* solution);
 		void  ResetConstraint(void);
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Node.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Node.cpp	(revision 16227)
+++ /issm/trunk-jpl/src/c/classes/Node.cpp	(revision 16228)
@@ -574,16 +574,4 @@
 	return indexing.clone;
 
-}
-/*}}}*/
-/*FUNCTION Node::UpdateSpcs {{{*/
-void   Node::UpdateSpcs(IssmDouble* ys){
-
-	int count = 0;
-	for(int i=0;i<this->indexing.gsize;i++){
-		if(this->indexing.s_set[i]){
-			this->indexing.svalues[i]=ys[this->indexing.sdoflist[count]];
-			count++;
-		}
-	}
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Node.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Node.h	(revision 16227)
+++ /issm/trunk-jpl/src/c/classes/Node.h	(revision 16228)
@@ -75,5 +75,4 @@
 		void  Activate(void);
 		void  Deactivate(void);
-		void  UpdateSpcs(IssmDouble* ys);
 		void  ReindexingDone(void);
 		bool  RequiresDofReindexing(void);
Index: /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp	(revision 16227)
+++ /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp	(revision 16228)
@@ -587,14 +587,2 @@
 }
 /*}}}*/
-/*FUNCTION ElementMatrix::SetDiag{{{*/
-void ElementMatrix::SetDiag(IssmDouble scalar){
-
-	int i;
-
-	if(this->nrows!=this->ncols)_error_("need square matrix in input!");
-
-	for(i=0;i<this->nrows;i++){
-		this->values[this->ncols*i+i]=scalar;
-	}
-}
-/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.h	(revision 16227)
+++ /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.h	(revision 16228)
@@ -65,5 +65,4 @@
 		void Transpose(void);
 		void Init(ElementMatrix* Ke);
-		void SetDiag(IssmDouble scalar);
 };
 #endif //#ifndef _ELEMENT_MATRIX_H_
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp	(revision 16227)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateDataSets.cpp	(revision 16228)
@@ -182,7 +182,12 @@
 	}
 
-	/*Update Elements and Materials For Control methods*/
+	/*Update Elements and Materials For Inversions*/
 	#ifdef _HAVE_CONTROL_
 	UpdateElementsAndMaterialsControl(elements,materials,iomodel);
+	#endif
+
+	/*Update Elements and Materials For Dakota*/
+	#ifdef _HAVE_DAKOTA_
+	UpdateElementsAndMaterialsDakota(elements,materials,iomodel);
 	#endif
 
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Stressbalance/UpdateElementsStressbalance.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Stressbalance/UpdateElementsStressbalance.cpp	(revision 16227)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Stressbalance/UpdateElementsStressbalance.cpp	(revision 16228)
@@ -113,5 +113,5 @@
 		if(dakota_analysis)elements->InputDuplicate(VzEnum,QmuVzEnum);
 		if(isFS){
-			iomodel->FetchDataToInput(elements,PressureEnum);
+			iomodel->FetchDataToInput(elements,PressureEnum,0.);
 			if(dakota_analysis)elements->InputDuplicate(PressureEnum,QmuPressureEnum);
 		}
Index: /issm/trunk-jpl/src/c/shared/Exceptions/exceptions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Exceptions/exceptions.h	(revision 16227)
+++ /issm/trunk-jpl/src/c/shared/Exceptions/exceptions.h	(revision 16228)
@@ -92,6 +92,3 @@
 };
 /*}}}*/
-
-char* exprintf(const char* format,...);
-
 #endif
Index: sm/trunk-jpl/src/c/shared/Exceptions/exprintf.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Exceptions/exprintf.cpp	(revision 16227)
+++ 	(revision )
@@ -1,48 +1,0 @@
-/*!\file:  exprintf
- * \brief this is a modification of the sprintf function. 
- * Instead of returning an int, it will return the char* itself.
- * The advantage is to be able to do things like: 
- * ErrorException(exprintf("%s%i\n","test failed for id:",id));
- */ 
-
-#include <stdarg.h>
-#include <stdio.h>
-#include "../MemOps/MemOps.h"
-
-char* exprintf(const char* format,...){
-
-	/*returned string: */
-	char *buffer = NULL;
-	int   n,size = 100;
-
-	//variable list of arguments
-	va_list args;
-
-	while(true){
-
-		/*allocate buffer for given string size*/
-		buffer=xNew<char>(size);
-
-		/* Try to print in the allocated space. */
-		va_start(args, format);
-#ifndef WIN32
-		n=vsnprintf(buffer,size,format,args);
-#else
-		n=vsnprintf(buffer,size,format,args);
-#endif
-		va_end(args);
-
-		/* If that worked, return the string. */
-		if(n>-1 && n<size) break;
-
-		/* Else try again with more space. */
-		if(n>-1)   /* glibc 2.1 */
-		 size=n+1; /* precisely what is needed */
-		else       /* glibc 2.0 */
-		 size*=2;  /* twice the old size */
-
-		xDelete<char>(buffer);
-	}
-
-	return buffer;
-}
Index: sm/trunk-jpl/src/c/toolkits/petsc/patches/VecDuplicatePatch.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/patches/VecDuplicatePatch.cpp	(revision 16227)
+++ 	(revision )
@@ -1,21 +1,0 @@
-/*!\file:  VecDuplicatePatch.cpp
- * \brief VecDuplicate + VecCopy wrapped together
- */ 
-
-#ifdef HAVE_CONFIG_H
-	#include <config.h>
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-/*Petsc includes: */
-#include <petscmat.h>
-#include <petscvec.h>
-#include <petscksp.h>
-
-void VecDuplicatePatch(Vec* output, Vec input){
-
-	VecDuplicate(input,output);
-	VecCopy(input,*output);
-
-}
Index: /issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h	(revision 16227)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h	(revision 16228)
@@ -34,5 +34,4 @@
 void MatMultPatch(Mat A,Vec X, Vec AX,ISSM_MPI_Comm comm);
 void MatToSerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm);
-void VecDuplicatePatch(Vec* output, Vec input);
 Vec  SerialToVec(double* vector,int vector_size);
 InsertMode ISSMToPetscInsertMode(InsMode mode);
