Index: /issm/trunk/src/c/objects/Loads/Numericalflux.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Numericalflux.cpp	(revision 3883)
+++ /issm/trunk/src/c/objects/Loads/Numericalflux.cpp	(revision 3884)
@@ -63,4 +63,5 @@
 	int  num_nodes;
 	int  num_elems;
+	bool   isinternal;
 
 	/*numericalflux constructor data: */
@@ -73,9 +74,9 @@
 	numericalflux_mparid=iomodel->numberofelements+1; //matlab indexing
 
-	/*Get left and right elements*/
-	e1=(int)iomodel->edges[4*i+2]; //edges are [node1 node2 elem1 elem2]
-	e2=(int)iomodel->edges[4*i+3]; //edges are [node1 node2 elem1 elem2]
-	if (isnan(e2)){
+	/*First, see wether this is an internal or boundary edge (if e2=NaN)*/
+	if (isnan((double)iomodel->edges[4*i+3])){ //edges are [node1 node2 elem1 elem2]
 		/* Boundary edge, only one element */
+		e1=(int)iomodel->edges[4*i+2];
+		e2=(int)UNDEF;
 		num_elems=1;
 		num_nodes=2;
@@ -85,4 +86,6 @@
 	else{
 		/* internal edge: connected to 2 elements */
+		e1=(int)iomodel->edges[4*i+2];
+		e2=(int)iomodel->edges[4*i+3];
 		num_elems=2;
 		num_nodes=4;
@@ -103,6 +106,6 @@
 		pos1=pos2=UNDEF;
 		for(j=0;j<3;j++){
-			if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1;
-			if (iomodel->elements[3*(int)e2+j]==i1) pos2=j+1;
+			if (iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
+			if (iomodel->elements[3*(e2-1)+j]==i1) pos2=j+1;
 		}
 		ISSMASSERT(pos1!=UNDEF && pos2!=UNDEF);
@@ -110,8 +113,8 @@
 		/*3: We have the id of the elements and the position of the vertices in the index
 		 * we can compute their dofs!*/
-		numericalflux_node_ids[0]=3*(int)e1+pos1;       //ex: 1 2 3
-		numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1; //ex: 2 3 1
-		numericalflux_node_ids[2]=3*(int)e2+pos2;           //ex: 1 2 3
-		numericalflux_node_ids[3]=3*(int)e2+((pos2+1)%3)+1; //ex: 3 1 2
+		numericalflux_node_ids[0]=3*(e1-1)+pos1;       //ex: 1 2 3
+		numericalflux_node_ids[1]=3*(e1-1)+(pos1%3)+1; //ex: 2 3 1
+		numericalflux_node_ids[2]=3*(e2-1)+pos2;           //ex: 1 2 3
+		numericalflux_node_ids[3]=3*(e2-1)+((pos2+1)%3)+1; //ex: 3 1 2
 	}
 	else{
@@ -120,5 +123,5 @@
 		pos1==UNDEF;
 		for(j=0;j<3;j++){
-			if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1;
+			if (iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
 		}
 		ISSMASSERT(pos1!=UNDEF);
@@ -126,7 +129,6 @@
 		/*3: We have the id of the elements and the position of the vertices in the index
 		 * we can compute their dofs!*/
-		numericalflux_node_ids[0]=3*(int)e1+pos1;
-		numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1;
-
+		numericalflux_node_ids[0]=3*(e1-1)+pos1;
+		numericalflux_node_ids[1]=3*(e1-1)+(pos1%3)+1;
 	}