Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 23879)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 23880)
@@ -2600,13 +2600,8 @@
 	IssmDouble mass_flux=0.;
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble normal[2];
-	IssmDouble length,rho_ice;
-	IssmDouble h1,h2;
-	IssmDouble vx1,vx2,vy1,vy2;
-	GaussTria* gauss_1=NULL;
-	GaussTria* gauss_2=NULL;
+	IssmDouble vx1,vx2,vy1,vy2,h1,h2;
 
 	/*Get material parameters :*/
-	rho_ice=FindParam(MaterialsRhoIceEnum);
+	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
 
 	/*First off, check that this segment belongs to this element: */
@@ -2617,16 +2612,17 @@
 
 	/*get area coordinates of 0 and 1 locations: */
-	gauss_1=new GaussTria();
+	GaussTria* gauss_1=new GaussTria();
 	gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
-	gauss_2=new GaussTria();
+	GaussTria* gauss_2=new GaussTria();
 	gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
 
-	normal[0]=cos(atan2(x1-x2,y2-y1));
-	normal[1]=sin(atan2(x1-x2,y2-y1));
-
-	length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
-
+	/*Get segment length and normal (needs to be properly oriented)*/
+	IssmDouble nx=cos(atan2(x1-x2,y2-y1));
+	IssmDouble ny=sin(atan2(x1-x2,y2-y1));
+	IssmDouble length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
+
+	/*Get velocity and thickness*/
+	this->parameters->FindParam(&domaintype,DomainTypeEnum);
 	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
-	this->parameters->FindParam(&domaintype,DomainTypeEnum);
 	Input* vx_input=NULL;
 	Input* vy_input=NULL;
@@ -2648,6 +2644,6 @@
 
 	mass_flux= rho_ice*length*(  
-				(ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
-				(ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
+				(1./3.*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*nx+
+				(1./3.*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*ny
 				);
 
@@ -2659,67 +2655,5 @@
 /*}}}*/
 IssmDouble Tria::MassFlux(IssmDouble* segment){/*{{{*/
-
-	int        domaintype;
-	IssmDouble mass_flux=0.;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble normal[2];
-	IssmDouble length,rho_ice;
-	IssmDouble x1,y1,x2,y2,h1,h2;
-	IssmDouble vx1,vx2,vy1,vy2;
-	GaussTria* gauss_1=NULL;
-	GaussTria* gauss_2=NULL;
-
-	/*Get material parameters :*/
-	rho_ice=FindParam(MaterialsRhoIceEnum);
-
-	/*First off, check that this segment belongs to this element: */
-	if (reCast<int>(*(segment+4))!=this->id)_error_("error message: segment with id " << reCast<int>(*(segment+4)) << " does not belong to element with id:" << this->id);
-
-	/*Recover segment node locations: */
-	x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
-
-	/*Get xyz list: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*get area coordinates of 0 and 1 locations: */
-	gauss_1=new GaussTria();
-	gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
-	gauss_2=new GaussTria();
-	gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
-
-	normal[0]=cos(atan2(x1-x2,y2-y1));
-	normal[1]=sin(atan2(x1-x2,y2-y1));
-
-	length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
-
-	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
-	this->parameters->FindParam(&domaintype,DomainTypeEnum);
-	Input* vx_input=NULL;
-	Input* vy_input=NULL;
-	if(domaintype==Domain2DhorizontalEnum){
-		vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
-		vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
-	}
-	else{
-		vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
-		vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
-	}
-
-	thickness_input->GetInputValue(&h1, gauss_1);
-	thickness_input->GetInputValue(&h2, gauss_2);
-	vx_input->GetInputValue(&vx1,gauss_1);
-	vx_input->GetInputValue(&vx2,gauss_2);
-	vy_input->GetInputValue(&vy1,gauss_1);
-	vy_input->GetInputValue(&vy2,gauss_2);
-
-	mass_flux= rho_ice*length*(  
-				(ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
-				(ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
-				);
-
-	/*clean up and return:*/
-	delete gauss_1;
-	delete gauss_2;
-	return mass_flux;
+	return this->MassFlux(segment[0],segment[1],segment[2],segment[3],reCast<int>(segment[4]));
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/shared/Numerics/constants.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/constants.h	(revision 23879)
+++ /issm/trunk-jpl/src/c/shared/Numerics/constants.h	(revision 23880)
@@ -7,5 +7,4 @@
 
 #define UNDEF -9999
-#define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
 #define SQRT2 1.414213562373095048801688724209698078569671875376948073176679738
 #define SQRT3 1.732050807568877293527446341505872366942805253810380628055806979