Index: /issm/trunk/src/c/objects/Element.h
===================================================================
--- /issm/trunk/src/c/objects/Element.h	(revision 3564)
+++ /issm/trunk/src/c/objects/Element.h	(revision 3565)
@@ -42,5 +42,4 @@
 		int            Enum();
 
-		
 };
 #endif
Index: /issm/trunk/src/c/objects/Numericalflux.cpp
===================================================================
--- /issm/trunk/src/c/objects/Numericalflux.cpp	(revision 3564)
+++ /issm/trunk/src/c/objects/Numericalflux.cpp	(revision 3565)
@@ -154,244 +154,18 @@
 void  Numericalflux::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
 
-	/*No stiffness loads applied, do nothing: */
-	return;
-
-}
-/*}}}*/
-/*FUNCTION Numericalflux::CreatePVector {{{1*/
-void  Numericalflux::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
-
-	/*No stiffness loads applied, do nothing: */
-	return;
-}
-/*}}}*/
-/*FUNCTION Numericalflux::DeepEcho {{{1*/
-void Numericalflux::DeepEcho(void){
-
-	int i;
-	
-	printf("Numericalflux:\n");
-	printf("   type: %s\n",type);
-	printf("   id: %i\n",id);
-	
-	printf("   element_id=%i\n",element_id);
-	printf("   element_offset=%i\n",element_offset);
-	if(element)element->Echo();
 	if (strcmp(type,"internal")==0){
-		printf("   node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]);
-		printf("   node_offsets=[%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]);
-		for(i=0;i<4;i++){
-			if(nodes[i])nodes[i]->Echo();
-		}
-	}
-	else{
-		printf("   node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]);
-		printf("   node_offsets=[%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]);
-		for(i=0;i<2;i++){
-			if(nodes[i])nodes[i]->Echo();
-		}
-	}
-	return;
-}		
-/*}}}*/
-/*FUNCTION Numericalflux::DistributeNumDofs {{{1*/
-void  Numericalflux::DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type){
-	return;
-}
-/*}}}*/
-/*FUNCTION Numericalflux::Echo {{{1*/
-void Numericalflux::Echo(void){
-
-	printf("Numericalflux:\n");
-	printf("   type: %s\n",type);
-	printf("   id: %i\n",id);
-
-	printf("   element_id=%i\n",element_id);
-	printf("   element_offset=%i]\n",element_offset);
-
-	if (strcmp(type,"internal")==0){
-		printf("   node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]);
-		printf("   node_offsets=[%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]);
-	}
-	else{
-		printf("   node_ids=[%i,%i]\n",node_ids[0],node_ids[1]);
-		printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1]);
-	}
-
-	return;
-}
-/*}}}*/
-/*FUNCTION Numericalflux::Enum {{{1*/
-int Numericalflux::Enum(void){
-
-	return NumericalfluxEnum();
-
-}
-/*}}}*/
-/*FUNCTION Numericalflux::GetB {{{1*/
-void Numericalflux::GetB(double* B, double gauss_coord){
-
-	const int numgrids=4;
-	double l1l4[numgrids];
-
-	/*Get nodal functions: */
-	GetNodalFunctions(&l1l4[0],gauss_coord);
-
-	/*Build B: */
-	B[0] = +l1l4[0];
-	B[1] = +l1l4[1];
-	B[2] = -l1l4[0];
-	B[3] = -l1l4[1];
-}
-/*}}}*/
-/*FUNCTION Numericalflux::GetDofList{{{1*/
-
-void  Numericalflux::GetDofList(int* doflist,int* pnumberofdofspernode){
-
-	int i,j;
-	int doflist_per_node[MAXDOFSPERNODE];
-	int numberofdofspernode;
-	
-	if (strcmp(type,"internal")==0){
-		for(i=0;i<4;i++){
-			nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
-			for(j=0;j<numberofdofspernode;j++){
-				doflist[i*numberofdofspernode+j]=doflist_per_node[j];
-			}
-		}
+
+		CreateKMatrixInternal(Kgg,inputs,analysis_type,sub_analysis_type);
 	}
 	else if (strcmp(type,"boundary")==0){
-		for(i=0;i<2;i++){
-			nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
-			for(j=0;j<numberofdofspernode;j++){
-				doflist[i*numberofdofspernode+j]=doflist_per_node[j];
-			}
-		}
-	}
-	else ISSMERROR(exprintf("type %s not supported yet",type));
-
-	/*Assign output pointers:*/
-	*pnumberofdofspernode=numberofdofspernode;
-}
-/*}}}*/
-/*FUNCTION Numericalflux::GetId {{{1*/
-int    Numericalflux::GetId(void){
-	return id;
-}
-/*}}}*/
-/*FUNCTION Numericalflux::GetJacobianDeterminant{{{1*/
-void Numericalflux::GetJacobianDeterminant(double* pJdet,double xyz_list[4][3], double gauss_coord){
-
-	double Jdet,length;
-
-	length=sqrt(pow(xyz_list[1][0] - xyz_list[0][0],2.0) + pow(xyz_list[1][1] - xyz_list[0][1],2.0));
-	Jdet=1.0/2.0*length;
-
-	if(Jdet<0){
-		ISSMERROR(" negative jacobian determinant!");
-	}
-
-	*pJdet=Jdet;
-}
-/*}}}*/
-/*FUNCTION Numericalflux::GetL {{{1*/
-void Numericalflux::GetL(double* L, double gauss_coord, int numdof){
-
-	const int numgrids=4;
-	double l1l4[numgrids];
-
-	/*Check numdof*/
-	ISSMASSERT(numdof==2 || numdof==4);
-
-	/*Get nodal functions: */
-	GetNodalFunctions(&l1l4[0],gauss_coord);
-
-	/*Luild L: */
-	L[0] = l1l4[0];
-	L[1] = l1l4[1];
-	if (numdof==4){
-		L[2] = l1l4[2];
-		L[3] = l1l4[3];
-	}
-}
-/*}}}*/
-/*FUNCTION Numericalflux::GetName {{{1*/
-char* Numericalflux::GetName(void){
-	return "numericalflux";
-}
-/*}}}*/
-/*FUNCTION Numericalflux::GetNodalFunctions{{{1*/
-void Numericalflux::GetNodalFunctions(double* l1l4, double gauss_coord){
-
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	l1l4[0]=-0.5*gauss_coord+0.5;
-	l1l4[1]=+0.5*gauss_coord+0.5;
-	l1l4[2]=-0.5*gauss_coord+0.5;
-	l1l4[3]=+0.5*gauss_coord+0.5;
-
-}
-/*}}}*/
-/*FUNCTION Numericalflux::GetNormal {{{1*/
-void Numericalflux:: GetNormal(double* normal,double xyz_list[4][3]){
-
-	/*Build unit outward pointing vector*/
-	const int numgrids=4;
-	double vector[2];
-	double norm;
-
-	vector[0]=xyz_list[1][0] - xyz_list[0][0];
-	vector[1]=xyz_list[1][1] - xyz_list[0][1];
-
-	norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
-
-	normal[0]= + vector[1]/norm;
-	normal[1]= - vector[0]/norm;
-}
-/*}}}*/
-/*FUNCTION Numericalflux::GetParameterValue {{{1*/
-void Numericalflux::GetParameterValue(double* pp, double* plist, double gauss_coord){
-
-	/*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian 
-	 * point specifie by gauss_l1l2l3: */
-
-	/*nodal functions: */
-	double l1l4[4];
-
-	/*output: */
-	double p;
-
-	GetNodalFunctions(&l1l4[0],gauss_coord);
-
-	p=l1l4[0]*plist[0]+l1l4[1]*plist[1];
-
-	/*Assign output pointers:*/
-	*pp=p;
-}
-/*}}}*/
-/*FUNCTION Numericalflux::MyRank {{{1*/
-int    Numericalflux::MyRank(void){ 
-	extern int my_rank;
-	return my_rank; 
-}
-/*}}}*/
-/*FUNCTION Numericalflux::PenaltyCreateKMatrix {{{1*/
-void  Numericalflux::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
-
-	if (strcmp(type,"internal")==0){
-
-		PenaltyCreateKMatrixInternal(Kgg,inputs,kmax,analysis_type,sub_analysis_type);
-	}
-	else if (strcmp(type,"boundary")==0){
-
-		PenaltyCreateKMatrixBoundary(Kgg,inputs,kmax,analysis_type,sub_analysis_type);
+
+		CreateKMatrixBoundary(Kgg,inputs,analysis_type,sub_analysis_type);
 	}
 	else ISSMERROR("type not supported yet");
-}
-/*}}}*/
-/*FUNCTION Numericalflux::PenaltyCreateKMatrixInternal {{{1*/
-void  Numericalflux::PenaltyCreateKMatrixInternal(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
-	double   start, finish;
-	start=MPI_Wtime();
+
+}
+/*}}}*/
+/*FUNCTION Numericalflux::CreateKMatrixInternal {{{1*/
+void  Numericalflux::CreateKMatrixInternal(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
 
 	/* local declarations */
@@ -494,8 +268,5 @@
 
 	/*Add Ke_gg to global matrix Kgg: */
-	start=MPI_Wtime();
 	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
-	finish=MPI_Wtime();
-	//printf("internal inserting K of edge %i matrix: %g\n",id,finish-start);
 
 	xfree((void**)&gauss_coords);
@@ -503,8 +274,6 @@
 }
 /*}}}*/
-/*FUNCTION Numericalflux::PenaltyCreateKMatrixBoundary {{{1*/
-void  Numericalflux::PenaltyCreateKMatrixBoundary(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
-	double   start, finish;
-	start=MPI_Wtime();
+/*FUNCTION Numericalflux::CreateKMatrixBoundary {{{1*/
+void  Numericalflux::CreateKMatrixBoundary(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
 
 	/* local declarations */
@@ -609,8 +378,5 @@
 
 	/*Add Ke_gg to global matrix Kgg: */
-	start=MPI_Wtime();
 	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
-	finish=MPI_Wtime();
-	//printf("boundary inserting K of edge %i matrix: %g\n",id,finish-start);
 
 	xfree((void**)&gauss_coords);
@@ -619,20 +385,21 @@
 }
 /*}}}*/
-/*FUNCTION Numericalflux::PenaltyCreatePVector{{{1*/
-void  Numericalflux::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
+/*FUNCTION Numericalflux::CreatePVector {{{1*/
+void  Numericalflux::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
 
 	if (strcmp(type,"internal")==0){
 
-		PenaltyCreatePVectorInternal(pg,inputs,kmax,analysis_type,sub_analysis_type);
+		CreatePVectorInternal(pg,inputs,analysis_type,sub_analysis_type);
 	}
 	else if (strcmp(type,"boundary")==0){
 
-		PenaltyCreatePVectorBoundary(pg,inputs,kmax,analysis_type,sub_analysis_type);
+		CreatePVectorBoundary(pg,inputs,analysis_type,sub_analysis_type);
 	}
 	else ISSMERROR("type not supported yet");
-}
-/*}}}*/
-/*FUNCTION Numericalflux::PenaltyCreatePVectorInternal{{{1*/
-void  Numericalflux::PenaltyCreatePVectorInternal(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
+
+}
+/*}}}*/
+/*FUNCTION Numericalflux::CreatePVectorInternal{{{1*/
+void  Numericalflux::CreatePVectorInternal(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
 
 	/*Nothing added to PVector*/
@@ -641,6 +408,6 @@
 }
 /*}}}*/
-/*FUNCTION Numericalflux::PenaltyCreatePVectorBoundary{{{1*/
-void  Numericalflux::PenaltyCreatePVectorBoundary(Vec pg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
+/*FUNCTION Numericalflux::CreatePVectorBoundary{{{1*/
+void  Numericalflux::CreatePVectorBoundary(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type){
 
 	/* local declarations */
@@ -745,4 +512,230 @@
 }
 /*}}}*/
+/*FUNCTION Numericalflux::DeepEcho {{{1*/
+void Numericalflux::DeepEcho(void){
+
+	int i;
+	
+	printf("Numericalflux:\n");
+	printf("   type: %s\n",type);
+	printf("   id: %i\n",id);
+	
+	printf("   element_id=%i\n",element_id);
+	printf("   element_offset=%i\n",element_offset);
+	if(element)element->Echo();
+	if (strcmp(type,"internal")==0){
+		printf("   node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]);
+		printf("   node_offsets=[%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]);
+		for(i=0;i<4;i++){
+			if(nodes[i])nodes[i]->Echo();
+		}
+	}
+	else{
+		printf("   node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]);
+		printf("   node_offsets=[%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]);
+		for(i=0;i<2;i++){
+			if(nodes[i])nodes[i]->Echo();
+		}
+	}
+	return;
+}		
+/*}}}*/
+/*FUNCTION Numericalflux::DistributeNumDofs {{{1*/
+void  Numericalflux::DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type){
+	return;
+}
+/*}}}*/
+/*FUNCTION Numericalflux::Echo {{{1*/
+void Numericalflux::Echo(void){
+
+	printf("Numericalflux:\n");
+	printf("   type: %s\n",type);
+	printf("   id: %i\n",id);
+
+	printf("   element_id=%i\n",element_id);
+	printf("   element_offset=%i]\n",element_offset);
+
+	if (strcmp(type,"internal")==0){
+		printf("   node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]);
+		printf("   node_offsets=[%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]);
+	}
+	else{
+		printf("   node_ids=[%i,%i]\n",node_ids[0],node_ids[1]);
+		printf("   node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1]);
+	}
+
+	return;
+}
+/*}}}*/
+/*FUNCTION Numericalflux::Enum {{{1*/
+int Numericalflux::Enum(void){
+
+	return NumericalfluxEnum();
+
+}
+/*}}}*/
+/*FUNCTION Numericalflux::GetB {{{1*/
+void Numericalflux::GetB(double* B, double gauss_coord){
+
+	const int numgrids=4;
+	double l1l4[numgrids];
+
+	/*Get nodal functions: */
+	GetNodalFunctions(&l1l4[0],gauss_coord);
+
+	/*Build B: */
+	B[0] = +l1l4[0];
+	B[1] = +l1l4[1];
+	B[2] = -l1l4[0];
+	B[3] = -l1l4[1];
+}
+/*}}}*/
+/*FUNCTION Numericalflux::GetDofList{{{1*/
+
+void  Numericalflux::GetDofList(int* doflist,int* pnumberofdofspernode){
+
+	int i,j;
+	int doflist_per_node[MAXDOFSPERNODE];
+	int numberofdofspernode;
+	
+	if (strcmp(type,"internal")==0){
+		for(i=0;i<4;i++){
+			nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
+			for(j=0;j<numberofdofspernode;j++){
+				doflist[i*numberofdofspernode+j]=doflist_per_node[j];
+			}
+		}
+	}
+	else if (strcmp(type,"boundary")==0){
+		for(i=0;i<2;i++){
+			nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
+			for(j=0;j<numberofdofspernode;j++){
+				doflist[i*numberofdofspernode+j]=doflist_per_node[j];
+			}
+		}
+	}
+	else ISSMERROR(exprintf("type %s not supported yet",type));
+
+	/*Assign output pointers:*/
+	*pnumberofdofspernode=numberofdofspernode;
+}
+/*}}}*/
+/*FUNCTION Numericalflux::GetId {{{1*/
+int    Numericalflux::GetId(void){
+	return id;
+}
+/*}}}*/
+/*FUNCTION Numericalflux::GetJacobianDeterminant{{{1*/
+void Numericalflux::GetJacobianDeterminant(double* pJdet,double xyz_list[4][3], double gauss_coord){
+
+	double Jdet,length;
+
+	length=sqrt(pow(xyz_list[1][0] - xyz_list[0][0],2.0) + pow(xyz_list[1][1] - xyz_list[0][1],2.0));
+	Jdet=1.0/2.0*length;
+
+	if(Jdet<0){
+		ISSMERROR(" negative jacobian determinant!");
+	}
+
+	*pJdet=Jdet;
+}
+/*}}}*/
+/*FUNCTION Numericalflux::GetL {{{1*/
+void Numericalflux::GetL(double* L, double gauss_coord, int numdof){
+
+	const int numgrids=4;
+	double l1l4[numgrids];
+
+	/*Check numdof*/
+	ISSMASSERT(numdof==2 || numdof==4);
+
+	/*Get nodal functions: */
+	GetNodalFunctions(&l1l4[0],gauss_coord);
+
+	/*Luild L: */
+	L[0] = l1l4[0];
+	L[1] = l1l4[1];
+	if (numdof==4){
+		L[2] = l1l4[2];
+		L[3] = l1l4[3];
+	}
+}
+/*}}}*/
+/*FUNCTION Numericalflux::GetName {{{1*/
+char* Numericalflux::GetName(void){
+	return "numericalflux";
+}
+/*}}}*/
+/*FUNCTION Numericalflux::GetNodalFunctions{{{1*/
+void Numericalflux::GetNodalFunctions(double* l1l4, double gauss_coord){
+
+	/*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+	l1l4[0]=-0.5*gauss_coord+0.5;
+	l1l4[1]=+0.5*gauss_coord+0.5;
+	l1l4[2]=-0.5*gauss_coord+0.5;
+	l1l4[3]=+0.5*gauss_coord+0.5;
+
+}
+/*}}}*/
+/*FUNCTION Numericalflux::GetNormal {{{1*/
+void Numericalflux:: GetNormal(double* normal,double xyz_list[4][3]){
+
+	/*Build unit outward pointing vector*/
+	const int numgrids=4;
+	double vector[2];
+	double norm;
+
+	vector[0]=xyz_list[1][0] - xyz_list[0][0];
+	vector[1]=xyz_list[1][1] - xyz_list[0][1];
+
+	norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
+
+	normal[0]= + vector[1]/norm;
+	normal[1]= - vector[0]/norm;
+}
+/*}}}*/
+/*FUNCTION Numericalflux::GetParameterValue {{{1*/
+void Numericalflux::GetParameterValue(double* pp, double* plist, double gauss_coord){
+
+	/*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian 
+	 * point specifie by gauss_l1l2l3: */
+
+	/*nodal functions: */
+	double l1l4[4];
+
+	/*output: */
+	double p;
+
+	GetNodalFunctions(&l1l4[0],gauss_coord);
+
+	p=l1l4[0]*plist[0]+l1l4[1]*plist[1];
+
+	/*Assign output pointers:*/
+	*pp=p;
+}
+/*}}}*/
+/*FUNCTION Numericalflux::MyRank {{{1*/
+int    Numericalflux::MyRank(void){ 
+	extern int my_rank;
+	return my_rank; 
+}
+/*}}}*/
+/*FUNCTION Numericalflux::PenaltyCreateKMatrix {{{1*/
+void  Numericalflux::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
+
+	/*No stiffness loads applied, do nothing: */
+	return;
+
+}
+/*}}}*/
+/*FUNCTION Numericalflux::PenaltyCreatePVector{{{1*/
+void  Numericalflux::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
+
+	/*No penalty loads applied, do nothing: */
+	return;
+
+}
+/*}}}*/
 /*FUNCTION Numericalflux::UpdateFromInputs {{{1*/
 void  Numericalflux::UpdateFromInputs(void* vinputs){
Index: /issm/trunk/src/c/objects/Numericalflux.h
===================================================================
--- /issm/trunk/src/c/objects/Numericalflux.h	(revision 3564)
+++ /issm/trunk/src/c/objects/Numericalflux.h	(revision 3565)
@@ -61,11 +61,12 @@
 		
 		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreateKMatrixInternal(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreateKMatrixBoundary(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type);
 		void  CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type);
+		void  CreatePVectorInternal(Vec pg,void* inputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVectorBoundary(Vec pg,void* inputs,int analysis_type,int sub_analysis_type);
 		void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
-		void  PenaltyCreateKMatrixInternal(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
-		void  PenaltyCreateKMatrixBoundary(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
 		void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
-		void  PenaltyCreatePVectorInternal(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
-		void  PenaltyCreatePVectorBoundary(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type);
+
 };