Index: /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp
===================================================================
--- /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp	(revision 3358)
+++ /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp	(revision 3359)
@@ -21,7 +21,11 @@
 
 
+	_printf_("   Configuring elements...\n");
 	elements->Configure(elements,loads,nodes,materials,parameters);
+	_printf_("   Configuring loads...\n");
 	loads->Configure(elements,loads,nodes,materials,parameters);
+	_printf_("   Configuring nodes...\n");
 	nodes->Configure(elements,loads,nodes,materials,parameters);
+	_printf_("   Configuring parameters...\n");
 	parameters->Configure(elements,loads, nodes, materials,parameters);
 
Index: /issm/trunk/src/c/DataSet/DataSet.cpp
===================================================================
--- /issm/trunk/src/c/DataSet/DataSet.cpp	(revision 3358)
+++ /issm/trunk/src/c/DataSet/DataSet.cpp	(revision 3359)
@@ -266,4 +266,10 @@
 			dataset->AddObject(icefront);
 		}
+		else if(enum_type==NumericalfluxEnum()){
+			Numericalflux* numericalflux=NULL;
+			numericalflux=new Numericalflux();
+			numericalflux->Demarshall(&marshalled_dataset);
+			dataset->AddObject(numericalflux);
+		}
 		else if(enum_type==RgbEnum()){
 			Rgb* rgb=NULL;
@@ -816,7 +822,5 @@
 		}
 		if(EnumIsLoad((*object)->Enum())){
-
 			load=(Load*)(*object);
-
 			load->Configure(elements,nodes,materials);
 		}
Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 3358)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 3359)
@@ -80,4 +80,5 @@
 int PenpairEnum(void){                      return          433; }
 int PengridEnum(void){                      return          434; }
+int NumericalfluxEnum(void){                return          435; }
 /*Materials: */
 int MaterialEnum(void){                     return          440; }
Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 3358)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 3359)
@@ -81,4 +81,5 @@
 int PenpairEnum(void);
 int PengridEnum(void);
+int NumericalfluxEnum(void);
 /*Materials: */
 int MaterialEnum(void);
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 3358)
+++ /issm/trunk/src/c/Makefile.am	(revision 3359)
@@ -84,4 +84,6 @@
 					./objects/Riftfront.cpp\
 					./objects/Riftfront.h\
+					./objects/Numericalflux.cpp\
+					./objects/Numericalflux.h\
 					./objects/Param.cpp\
 					./objects/Param.h\
@@ -468,4 +470,6 @@
 					./objects/Riftfront.cpp\
 					./objects/Riftfront.h\
+					./objects/Numericalflux.cpp\
+					./objects/Numericalflux.h\
 					./objects/Param.cpp\
 					./objects/Param.h\
Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 3358)
+++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 3359)
@@ -241,5 +241,6 @@
 	count++;
 	param= new Param(count,"numberofnodes",DOUBLE);
-	param->SetDouble(iomodel->numberofnodes);
+	if (iomodel->analysis_type==Prognostic2AnalysisEnum()) param->SetDouble(3*iomodel->numberofelements);
+	else param->SetDouble(iomodel->numberofnodes);
 	parameters->AddObject(param);
 
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp	(revision 3358)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp	(revision 3359)
@@ -95,5 +95,4 @@
 		}
 
-
 		element=(int)(*(iomodel->pressureload+segment_width*i+segment_width-2)-1); //element is in the penultimate column (grid1 grid2 ... elem fill)
 
Index: /issm/trunk/src/c/ModelProcessorx/IoModel.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/IoModel.cpp	(revision 3358)
+++ /issm/trunk/src/c/ModelProcessorx/IoModel.cpp	(revision 3359)
@@ -95,4 +95,6 @@
 	iomodel-> spctemperature=NULL;
 	iomodel-> spcthickness=NULL;
+	iomodel->numberofedges=0;
+	iomodel->edges=NULL;
 	
 	/*!materials: */
@@ -247,4 +249,5 @@
 	xfree((void**)&iomodel->spcthickness);
 	xfree((void**)&iomodel->spctemperature);
+	xfree((void**)&iomodel->edges);
 	xfree((void**)&iomodel->geothermalflux);
 	xfree((void**)&iomodel->melting);
Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/IoModel.h	(revision 3358)
+++ /issm/trunk/src/c/ModelProcessorx/IoModel.h	(revision 3359)
@@ -95,4 +95,6 @@
 	double* spcthickness;
 	double* geothermalflux;
+	int     numberofedges;
+	double* edges;
 	
 	/*materials: */
Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateConstraintsPrognostic2.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateConstraintsPrognostic2.cpp	(revision 3358)
+++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateConstraintsPrognostic2.cpp	(revision 3359)
@@ -2,5 +2,4 @@
  * CreateConstraintsPrognostic2.c:
  */
-
 
 #include "../../DataSet/DataSet.h"
@@ -11,64 +10,10 @@
 #include "../IoModel.h"
 
-
 void	CreateConstraintsPrognostic2(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
-
-
-	int i;
-	int count;
 	
 	DataSet* constraints = NULL;
-	Spc*    spc  = NULL;
-
-	/*spc intermediary data: */
-	int spc_sid;
-	int spc_node;
-	int spc_dof;
-	double spc_value;
-	
-	double* spcthickness=NULL;
 	
 	/*Create constraints: */
 	constraints = new DataSet(ConstraintsEnum());
-
-	/*Fetch data: */
-	IoModelFetchData(&spcthickness,NULL,NULL,iomodel_handle,"spcthickness");
-
-	count=0;
-
-	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
-	for (i=0;i<iomodel->numberofnodes;i++){
-	#ifdef _PARALLEL_
-	/*keep only this partition's nodes:*/
-	if((iomodel->my_grids[i]==1)){
-	#endif
-
-		if ((int)spcthickness[2*i]){
-	
-			/*This grid needs to be spc'd: */
-
-			spc_sid=count;
-			spc_node=i+1;
-			spc_dof=1; //we enforce first translation degree of freedom, for temperature
-			spc_value=*(spcthickness+2*i+1);
-
-			spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
-			constraints->AddObject(spc);
-			count++;
-		}
-
-	#ifdef _PARALLEL_
-	} //if((my_grids[i]==1))
-	#endif
-	}
-
-	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
-	 * datasets, it will not be redone: */
-	constraints->Presort();
-
-	/*Free data: */
-	xfree((void**)&spcthickness);
-	
-	cleanup_and_return:
 	
 	/*Assign output pointer: */
Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp	(revision 3358)
+++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp	(revision 3359)
@@ -246,82 +246,4 @@
 	else{ //	if (strcmp(meshtype,"2d")==0)
 		ISSMERROR(exprintf("not implemented yet"));
-
-		/*Fetch data needed: */
-		IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
-		IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
-		IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface");
-		IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
-		IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf");
-		IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed");
-		IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface");
-		IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
-	
-		for (i=0;i<iomodel->numberofelements;i++){
-		#ifdef _PARALLEL_
-		/*We are using our element partition to decide which elements will be created on this node: */
-		if(my_rank==epart[i]){
-		#endif
-
-			
-			/*name and id: */
-			penta_id=i+1; //matlab indexing.
-			penta_mid=-1; 
-			penta_mparid=-1; //no need for materials
-			penta_numparid=1;
-
-			/*vertices,thickness,surface,bed and drag: */
-			for(j=0;j<6;j++){
-				penta_g[j]=(int)*(iomodel->elements+elements_width*i+j);
-				penta_h[j]=*(iomodel->thickness+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
-				penta_s[j]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
-				penta_b[j]=*(iomodel->bed+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
-			}
-
-			/*diverse: */
-			penta_shelf=(int)*(iomodel->elementoniceshelf+i);
-			penta_onbed=(int)*(iomodel->elementonbed+i);
-			penta_onsurface=(int)*(iomodel->elementonsurface+i);
-			penta_collapse=1;
-			penta_onwater=(bool)*(iomodel->elementonwater+i);
-	
-
-			/*Create Penta using its constructor:*/
-			penta= new Penta( penta_id,penta_mid,penta_mparid,penta_numparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
-					penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,
-					penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,
-					penta_thermal_steadystate,penta_onwater); 
-
-			/*Add penta element to elements dataset: */
-			elements->AddObject(penta);
-	
-			#ifdef _PARALLEL_
-			/*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 
-			 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
-			 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 
-			 will hold which grids belong to this partition*/
-			my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
-			my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
-			my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
-			my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;
-			my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;
-			my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;
-			#endif
-
-		#ifdef _PARALLEL_
-		}//if(my_rank==epart[i])
-		#endif
-
-		}//for (i=0;i<numberofelements;i++)
-
-		/*Free data: */
-		xfree((void**)&iomodel->elements);
-		xfree((void**)&iomodel->thickness);
-		xfree((void**)&iomodel->surface);
-		xfree((void**)&iomodel->bed);
-		xfree((void**)&iomodel->elementoniceshelf);
-		xfree((void**)&iomodel->elementonbed);
-		xfree((void**)&iomodel->elementonsurface);
-		xfree((void**)&iomodel->elementonwater);
-
 	} //if (strcmp(meshtype,"2d")==0)
 
@@ -373,35 +295,4 @@
 	materials->AddObject(matpar);
 
-	/*Partition penalties in 3d: */
-	if(strcmp(iomodel->meshtype,"3d")==0){
-	
-		/*Get penalties: */
-		IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");
-
-		if(iomodel->numpenalties){
-
-			iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int));
-			#ifdef _SERIAL_
-			for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1;
-			#else
-			for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1;
-
-			for(i=0;i<iomodel->numpenalties;i++){
-				first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1);
-				if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition  grids
-					/*All grids that are being penalised belong to this node's internal grid partition.:*/
-					iomodel->penaltypartitioning[i]=1;
-				}
-				if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
-					iomodel->penaltypartitioning[i]=0;
-				}
-			}
-			#endif
-		}
-
-		/*Free penalties: */
-		xfree((void**)&iomodel->penalties);
-	}
-
 	/*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. 
 	 We can therefore determine  which grids are internal to this node's partition 
@@ -416,4 +307,5 @@
 		IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
 	}
+	IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
 	IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x");
 	IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y");
@@ -430,66 +322,50 @@
 	DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);
 
-	for (i=0;i<iomodel->numberofnodes;i++){
-	#ifdef _PARALLEL_
-	/*keep only this partition's nodes:*/
-	if((my_grids[i]==1)){
-	#endif
-
-		node_id=i+1; //matlab indexing
-			
-		
-		
-		#ifdef _PARALLEL_
-		if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
-			node_partitionborder=1;
-		}
-		else{
+	/*Build Nodes dataset -> 3 for each element*/
+	for (i=0;i<iomodel->numberofelements;i++){
+		for (j=0;j<3;j++){
+
+			#ifdef _PARALLEL_
+			/*!All elements have been partitioned above, only create elements for this CPU: */
+			if(my_rank==epart[i]){ 
+			#endif
+
+			//Get id of the vertex on which the current node is located
+			k=(int)*(iomodel->elements+elements_width*i+j)-1; //(Matlab to C indexing)
+
+			//Get node properties
+			node_id=i*3+j+1;
 			node_partitionborder=0;
-		}
-		#else
-			node_partitionborder=0;
-		#endif
-
-		node_x[0]=iomodel->x[i];
-		node_x[1]=iomodel->y[i];
-		node_x[2]=iomodel->z[i];
-		node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
-		
-		node_onbed=(int)iomodel->gridonbed[i];
-		node_onsurface=(int)iomodel->gridonsurface[i];	
-		node_onshelf=(int)iomodel->gridoniceshelf[i];	
-		node_onsheet=(int)iomodel->gridonicesheet[i];	
-
-		if (strcmp(iomodel->meshtype,"3d")==0){
-			if (isnan(iomodel->uppernodes[i])){
-				node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
+			node_x[0]=iomodel->x[k];
+			node_x[1]=iomodel->y[k];
+			node_x[2]=iomodel->z[k];
+			node_sigma=(iomodel->z[k]-iomodel->bed[k])/(iomodel->thickness[k]);
+			node_onbed=(int)iomodel->gridonbed[k];
+			node_onsurface=(int)iomodel->gridonsurface[k];	
+			node_onshelf=(int)iomodel->gridoniceshelf[k];	
+			node_onsheet=(int)iomodel->gridonicesheet[k];	
+			if (strcmp(iomodel->meshtype,"3d")==0){
+				if (isnan(iomodel->uppernodes[k])){
+					node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
+				}
+				else{
+					node_upper_node_id=(int)iomodel->uppernodes[k];
+				}
 			}
 			else{
-				node_upper_node_id=(int)iomodel->uppernodes[i];
+				/*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
+				node_upper_node_id=node_id;
 			}
+
+			/*Create node using its constructor: */
+			node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
+
+			/*Add node to nodes dataset: */
+			nodes->AddObject(node);
+
+			#ifdef _PARALLEL_
+			}
+			#endif
 		}
-		else{
-			/*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
-			node_upper_node_id=node_id;
-		}
-
-		/*Create node using its constructor: */
-		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
-
-		/*set single point constraints.: */
-		if (strcmp(iomodel->meshtype,"3d")==0){
-			/*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
-			if (!iomodel->gridonbed[i]){
-				for(k=1;k<=node_numdofs;k++){
-					node->FreezeDof(k);
-				}
-			}
-		}
-		/*Add node to nodes dataset: */
-		nodes->AddObject(node);
-
-	#ifdef _PARALLEL_
-	} //if((my_grids[i]==1))
-	#endif
 	}
 
@@ -502,4 +378,5 @@
 	/*Clean fetched data: */
 	xfree((void**)&iomodel->deadgrids);
+	xfree((void**)&iomodel->elements);
 	xfree((void**)&iomodel->x);
 	xfree((void**)&iomodel->y);
@@ -513,5 +390,4 @@
 	xfree((void**)&iomodel->gridoniceshelf);
 	
-
 	/*Keep partitioning information into iomodel*/
 	iomodel->epart=epart;
@@ -532,4 +408,3 @@
 	*pnodes=nodes;
 	*pmaterials=materials;
-
 }
Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp	(revision 3358)
+++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp	(revision 3359)
@@ -1,5 +1,4 @@
 /*! \file CreateLoadsPrognostic2.c:
  */
-
 
 #include "../../DataSet/DataSet.h"
@@ -9,20 +8,113 @@
 #include "../../shared/shared.h"
 #include "../../include/macros.h"
+#include "../../include/typedefs.h"
 #include "../IoModel.h"
-
 
 void	CreateLoadsPrognostic2(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
-	DataSet*    loads    = NULL;
+	int i,j;
+	int i1,i2;
+	int pos1,pos2;
+	double e1,e2;
+
+	extern int my_rank;
+	extern int num_procs;
+
+	DataSet*       loads        = NULL;
+	Numericalflux* numericalflux= NULL;
+
+	/*numericalflux intermediary data: */
+	char numericalflux_type[NUMERICALFLUXSTRING];
+	int  numericalflux_id;
+	int  numericalflux_node_ids[MAX_NUMERICALFLUX_NODES];
+	int  numericalflux_elem_ids[MAX_NUMERICALFLUX_ELEMS];
 
 	/*Create loads: */
 	loads   = new DataSet(LoadsEnum());
-	
-	
+
+	/*Get edges and elements*/
+	IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges");
+	IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+
+	/*First load data:*/
+	for (i=0;i<iomodel->numberofedges;i++){
+
+		/*Get left and right elements*/
+		e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2]
+		e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2]
+
+		#ifdef _PARALLEL_
+		if (iomodel->epart[e1]!=my_rank){
+			/*This load does not belong to this cluster node, as it references an element 
+			 *that does not belong to this node's partition. Drop this 'i':*/
+			continue;
+		}
+		#endif
+
+		/*Create load*/
+		numericalflux_id=i+1; //Matlab indexing
+
+		if (!isnan(e2)){
+			strcpy(numericalflux_type,"internal");
+
+			numericalflux_elem_ids[0]=(int)e1+1;//id is in matlab index
+			numericalflux_elem_ids[1]=(int)e2+1;//id is in matlab index
+
+			/*Now, we must get the nodes of the 4 nodes located on the edge*/
+
+			/*1: Get vertices ids*/
+			i1=(int)iomodel->edges[4*i+0];
+			i2=(int)iomodel->edges[4*i+1];
+
+			/*2: Get the column where these ids are located in the index*/
+			pos1=pos2=UNDEF;
+			for(j=0;j<3;j++){
+				if (iomodel->elements[3*(int)e1+j]==i1) pos1=j;
+				if (iomodel->elements[3*(int)e2+j]==i1) pos2=j;
+			}
+			ISSMASSERT(pos1!=UNDEF && pos2!=UNDEF);
+
+			/*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+1;      //id is in matlab index
+			numericalflux_node_ids[1]=3*(int)e1+(pos1+1)%3+1;
+			numericalflux_node_ids[2]=3*(int)e2+pos2+1;
+			numericalflux_node_ids[3]=3*(int)e2+(pos2-1)%3+1;
+		}
+		else{
+			strcpy(numericalflux_type,"boundary");
+
+			numericalflux_elem_ids[0]=(int)e1+1;
+
+			/*1: Get vertices ids*/
+			i1=(int)iomodel->edges[4*i+0];
+			i2=(int)iomodel->edges[4*i+1];
+
+			/*2: Get the column where these ids are located in the index*/
+			pos1==UNDEF;
+			for(j=0;j<3;j++){
+				if (iomodel->elements[3*(int)e1+j]==i1) pos1=j;
+			}
+			ISSMASSERT(pos1!=UNDEF);
+
+			/*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+1; //id is in matlab index
+			numericalflux_node_ids[1]=3*(int)e1+(pos1+1)%3+1;
+		}
+
+		numericalflux = new Numericalflux(numericalflux_id,numericalflux_type,numericalflux_node_ids,numericalflux_elem_ids);
+
+		loads->AddObject(numericalflux);
+
+	}
+
+	/*Free data: */
+	xfree((void**)&iomodel->edges);
+	xfree((void**)&iomodel->elements);
+
 	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
 	 * datasets, it will not be redone: */
 	loads->Presort();
-
-	cleanup_and_return:
 
 	/*Assign output pointer: */
@@ -30,4 +122,2 @@
 
 }
-
-
Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp	(revision 3358)
+++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp	(revision 3359)
@@ -9,4 +9,5 @@
 #include "../../objects/objects.h"
 #include "../../shared/shared.h"
+#include "../../include/macros.h"
 #include "../IoModel.h"
 
@@ -16,18 +17,19 @@
 	DataSet* parameters=NULL;
 	int      count;
-	int      i;
+	int      i,j;
 	int      dim;
+	double*  elements;
+	int*     part;
 
 	double* vx=NULL;
+	double* vx_g=NULL;
 	double* vy=NULL;
-	double* vz=NULL;
-	double* u_g=NULL;
-	double* pressure=NULL;
-	double* temperature=NULL;
+	double* vy_g=NULL;
 	double* thickness=NULL;
-	double* surface=NULL;
-	double* bed=NULL;
+	double* h_g=NULL;
 	double* accumulation=NULL;
+	double* a_g=NULL;
 	double* melting=NULL;
+	double* m_g=NULL;
 
 	/*recover parameters : */
@@ -36,118 +38,72 @@
 	count=parameters->Size();
 
-	/*Get vx and vy: */
+	/*Create transformatiob vector for DG inputs*/
+	IoModelFetchData(&elements,NULL,NULL,iomodel_handle,"elements");
+	part=(int*)xcalloc(iomodel->numberofelements*3,sizeof(int));
+	for(i=0;i<iomodel->numberofelements;i++){
+		for(j=0;j<3;j++){
+			part[3*i+j]=(int)elements[3*i+j]-1; //Matlab to C indexing
+			ISSMASSERT(part[3*i+j]<iomodel->numberofnodes);
+		}
+	}
+	xfree((void**)&elements);
+
+	/*Get vx: */
 	IoModelFetchData(&vx,NULL,NULL,iomodel_handle,"vx");
-	IoModelFetchData(&vy,NULL,NULL,iomodel_handle,"vy");
-	IoModelFetchData(&vz,NULL,NULL,iomodel_handle,"vz");
-
-	u_g=(double*)xcalloc(iomodel->numberofnodes*3,sizeof(double));
-
-	if(vx) for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+0]=vx[i]/iomodel->yts;
-	if(vy) for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+1]=vy[i]/iomodel->yts;
-	if(vz) for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+2]=vz[i]/iomodel->yts;
+	vx_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
+	if(vx) for(i=0;i<3*iomodel->numberofelements;i++) vx_g[i]=vx[part[i]]/iomodel->yts;
 
 	count++;
-	param= new Param(count,"u_g",DOUBLEVEC);
-	param->SetDoubleVec(u_g,3*iomodel->numberofnodes,3);
+	param= new Param(count,"vx_g",DOUBLEVEC);
+	param->SetDoubleVec(vx_g,3*iomodel->numberofelements,1);
 	parameters->AddObject(param);
+	xfree((void**)&vx);
 
+	/*Get vy: */
+	IoModelFetchData(&vy,NULL,NULL,iomodel_handle,"vy");
+	vy_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
+	if(vy) for(i=0;i<3*iomodel->numberofelements;i++) vy_g[i]=vy[part[i]]/iomodel->yts;
 
-	xfree((void**)&vx);
+	count++;
+	param= new Param(count,"vy_g",DOUBLEVEC);
+	param->SetDoubleVec(vy_g,3*iomodel->numberofelements,1);
+	parameters->AddObject(param);
 	xfree((void**)&vy);
-	xfree((void**)&vz);
-	xfree((void**)&u_g);
-
-	/*Get pressure if 3d iomodel: */
-	parameters->FindParam(&dim,"dim");
-	if (dim==3){ 
-		IoModelFetchData(&pressure,NULL,NULL,iomodel_handle,"pressure");
-		
-		count++;
-		param= new Param(count,"p_g",DOUBLEVEC);
-		if(pressure) param->SetDoubleVec(pressure,iomodel->numberofnodes,1);
-		else param->SetDoubleVec(pressure,0,0);
-		parameters->AddObject(param);
-
-		/*Free pressure: */
-		xfree((void**)&pressure);
-	}
-
-	/*Get temperature if 3d iomodel: */
-	parameters->FindParam(&dim,"dim");
-	if (dim==3){ 
-		IoModelFetchData(&temperature,NULL,NULL,iomodel_handle,"temperature");
-		
-		count++;
-		param= new Param(count,"t_g",DOUBLEVEC);
-		if(temperature) param->SetDoubleVec(temperature,iomodel->numberofnodes,1);
-		else param->SetDoubleVec(temperature,0,0);
-		parameters->AddObject(param);
-
-		/*Free temperature: */
-		xfree((void**)&temperature);
-	}
 
 	/*Get thickness: */
 	IoModelFetchData(&thickness,NULL,NULL,iomodel_handle,"thickness");
-	
+	h_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
+	if(thickness) for(i=0;i<3*iomodel->numberofelements;i++) h_g[i]=thickness[part[i]]/iomodel->yts;
+
 	count++;
 	param= new Param(count,"h_g",DOUBLEVEC);
-	if(thickness) param->SetDoubleVec(thickness,iomodel->numberofnodes,1);
-	else param->SetDoubleVec(thickness,0,0);
+	param->SetDoubleVec(h_g,3*iomodel->numberofelements,1);
 	parameters->AddObject(param);
-
-	/*Free thickness: */
 	xfree((void**)&thickness);
 
-	/*Get surface: */
-	IoModelFetchData(&surface,NULL,NULL,iomodel_handle,"surface");
-	
+	/*Get accumulation: */
+	IoModelFetchData(&accumulation,NULL,NULL,iomodel_handle,"accumulation");
+	a_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
+	if(accumulation) for(i=0;i<3*iomodel->numberofelements;i++) a_g[i]=accumulation[part[i]]/iomodel->yts;
+
 	count++;
-	param= new Param(count,"s_g",DOUBLEVEC);
-	if(surface) param->SetDoubleVec(surface,iomodel->numberofnodes,1);
-	else param->SetDoubleVec(surface,0,0);
+	param= new Param(count,"a_g",DOUBLEVEC);
+	param->SetDoubleVec(a_g,3*iomodel->numberofelements,1);
 	parameters->AddObject(param);
-
-	/*Free surface: */
-	xfree((void**)&surface);
-
-	/*Get bed: */
-	IoModelFetchData(&bed,NULL,NULL,iomodel_handle,"bed");
-	
-	count++;
-	param= new Param(count,"b_g",DOUBLEVEC);
-	if(bed) param->SetDoubleVec(bed,iomodel->numberofnodes,1);
-	else param->SetDoubleVec(bed,0,0);
-	parameters->AddObject(param);
-
-	/*Free bed: */
-	xfree((void**)&bed);
+	xfree((void**)&accumulation);
 
 	/*Get melting: */
 	IoModelFetchData(&melting,NULL,NULL,iomodel_handle,"melting");
-	if(melting) for(i=0;i<iomodel->numberofnodes;i++)melting[i]=melting[i]/iomodel->yts;
-	
+	m_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double));
+	if(melting) for(i=0;i<3*iomodel->numberofelements;i++) m_g[i]=melting[part[i]]/iomodel->yts;
+
 	count++;
 	param= new Param(count,"m_g",DOUBLEVEC);
-	if(melting) param->SetDoubleVec(melting,iomodel->numberofnodes,1);
-	else param->SetDoubleVec(melting,0,1);
+	param->SetDoubleVec(m_g,3*iomodel->numberofelements,1);
 	parameters->AddObject(param);
-
-	/*Free melting: */
 	xfree((void**)&melting);
 
-	/*Get accumulation: */
-	IoModelFetchData(&accumulation,NULL,NULL,iomodel_handle,"accumulation");
-	if(accumulation) for(i=0;i<iomodel->numberofnodes;i++)accumulation[i]=accumulation[i]/iomodel->yts;
-	
-	count++;
-	param= new Param(count,"a_g",DOUBLEVEC);
-	if(accumulation) param->SetDoubleVec(accumulation,iomodel->numberofnodes,1);
-	else param->SetDoubleVec(accumulation,0,0);
-	parameters->AddObject(param);
-
-	/*Free accumulation: */
-	xfree((void**)&accumulation);
-
+	/*Free partitioning vector*/
+	xfree((void**)&part);
 
 	/*Assign output pointer: */
Index: /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp	(revision 3358)
+++ /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp	(revision 3359)
@@ -14,14 +14,8 @@
 void	CreateConstraintsSlopeCompute(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
-
 	DataSet* constraints = NULL;
 
 	/*Create constraints: */
 	constraints = new DataSet(ConstraintsEnum());
-
-	/*Now, is the flag isstokes on? otherwise, do nothing: */
-	if (!iomodel->isstokes)goto cleanup_and_return;
-	
-	cleanup_and_return:	
 	
 	/*Assign output pointer: */
Index: /issm/trunk/src/c/include/macros.h
===================================================================
--- /issm/trunk/src/c/include/macros.h	(revision 3358)
+++ /issm/trunk/src/c/include/macros.h	(revision 3359)
@@ -20,5 +20,5 @@
 /*Assertion macro*/
 #define ISSMASSERT(statement)\
-  if (!(statement)) ISSMERROR(exprintf("Assertion %s failed, please report bug to an ISSM developer",#statement))
+  if (!(statement)) ISSMERROR(exprintf("Assertion \"%s\" failed, please report bug to an ISSM developer",#statement))
 
 /*The following macros hide the error exception handling in a matlab module. Just put 
Index: /issm/trunk/src/c/objects/Icefront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Icefront.cpp	(revision 3358)
+++ /issm/trunk/src/c/objects/Icefront.cpp	(revision 3359)
@@ -173,7 +173,5 @@
 	/*Link this load with its nodes: */
 	if (strcmp(type,"segment")==0){
-		
 		ResolvePointers((Object**)nodes,node_ids,node_offsets,2,nodesin);
-	
 	}
 	else{
Index: /issm/trunk/src/c/objects/Node.cpp
===================================================================
--- /issm/trunk/src/c/objects/Node.cpp	(revision 3358)
+++ /issm/trunk/src/c/objects/Node.cpp	(revision 3359)
@@ -600,4 +600,6 @@
 		}
 		else if ( 
+				(strcmp(field_name,"vx")==0) ||
+				(strcmp(field_name,"vy")==0) ||
 				(strcmp(field_name,"thickness")==0) ||
 				(strcmp(field_name,"temperature")==0) ||
Index: /issm/trunk/src/c/objects/Numericalflux.cpp
===================================================================
--- /issm/trunk/src/c/objects/Numericalflux.cpp	(revision 3359)
+++ /issm/trunk/src/c/objects/Numericalflux.cpp	(revision 3359)
@@ -0,0 +1,375 @@
+/*!\file Numericalflux.c
+ * \brief: implementation of the Numericalflux object
+ */
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "stdio.h"
+#include <string.h>
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../shared/shared.h"
+#include "../include/typedefs.h"
+#include "../include/macros.h"
+#include "./objects.h"
+
+/*Object constructors and destructor*/
+/*FUNCTION Numericalflux::Numericalflux(){{{1*/
+Numericalflux::Numericalflux(){
+	return;
+}
+/*}}}*/
+/*FUNCTION Numericalflux::Numericalflux(char numericalflux_type[NUMERICALFLUXSTRING],int numericalflux_fill...){{{1*/
+Numericalflux::Numericalflux(int numericalflux_id,char numericalflux_type[NUMERICALFLUXSTRING], int numericalflux_node_ids[MAX_NUMERICALFLUX_NODES],int numericalflux_element_ids[MAX_NUMERICALFLUX_ELEMS]){
+
+	int i;
+	
+	strcpy(type,numericalflux_type);
+	id=numericalflux_id;
+	
+	for(i=0;i<MAX_NUMERICALFLUX_ELEMS;i++){
+		element_ids[i]=numericalflux_element_ids[i];
+		element_offsets[i]=UNDEF;
+		elements[i]=NULL;
+	}
+
+	for(i=0;i<MAX_NUMERICALFLUX_NODES;i++){
+		node_ids[i]=numericalflux_node_ids[i];
+		node_offsets[i]=UNDEF;
+		nodes[i]=NULL;
+	}
+
+	return;
+}
+/*}}}*/
+/*FUNCTION Numericalflux::~Numericalflux(){{{1*/
+Numericalflux::~Numericalflux(){
+	return;
+}
+/*}}}*/
+/*FUNCTION Numericalflux::copy {{{1*/
+Object* Numericalflux::copy() {
+	return new Numericalflux(*this); 
+}
+/*}}}*/
+
+/*Object marshall*/
+/*FUNCTION Numericalflux::Demarshall {{{1*/
+void  Numericalflux::Demarshall(char** pmarshalled_dataset){
+
+	int i;
+	char* marshalled_dataset=NULL;
+
+	/*recover marshalled_dataset: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
+	 *object data (thanks to DataSet::Demarshall):*/
+
+	memcpy(&type,marshalled_dataset,sizeof(type));marshalled_dataset+=sizeof(type);
+	memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
+
+	memcpy(&element_ids,marshalled_dataset,sizeof(element_ids));marshalled_dataset+=sizeof(element_ids);
+	memcpy(&element_offsets,marshalled_dataset,sizeof(element_offsets));marshalled_dataset+=sizeof(element_offsets);
+	for(i=0;i<MAX_NUMERICALFLUX_ELEMS;i++)elements[i]=NULL;
+	
+	memcpy(&node_ids,marshalled_dataset,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
+	memcpy(&node_offsets,marshalled_dataset,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
+	for(i=0;i<MAX_NUMERICALFLUX_NODES;i++)nodes[i]=NULL;
+
+	/*return: */
+	*pmarshalled_dataset=marshalled_dataset;
+	return;
+}
+/*}}}*/
+/*FUNCTION Numericalflux::Marshall {{{1*/
+void  Numericalflux::Marshall(char** pmarshalled_dataset){
+
+	char* marshalled_dataset=NULL;
+	int   enum_type=0;
+
+	/*recover marshalled_dataset: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*get enum type of Numericalflux: */
+	enum_type=NumericalfluxEnum();
+	
+	/*marshall enum: */
+	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
+	
+	/*marshall Numericalflux data: */
+	memcpy(marshalled_dataset,&type,sizeof(type));marshalled_dataset+=sizeof(type);
+	memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
+	
+	memcpy(marshalled_dataset,&element_ids,sizeof(element_ids));marshalled_dataset+=sizeof(element_ids);
+	memcpy(marshalled_dataset,&element_offsets,sizeof(element_offsets));marshalled_dataset+=sizeof(element_offsets);
+
+	memcpy(marshalled_dataset,&node_ids,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
+	memcpy(marshalled_dataset,&node_offsets,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
+
+	*pmarshalled_dataset=marshalled_dataset;
+	return;
+}
+/*}}}*/
+/*FUNCTION Numericalflux::MarshallSize{{{1*/
+int   Numericalflux::MarshallSize(){
+
+	return sizeof(type)+
+		sizeof(id)+
+		sizeof(element_ids)+
+		sizeof(element_offsets)+
+		sizeof(node_ids)+
+		sizeof(node_offsets)+
+		sizeof(int); //sizeof(int) for enum type
+}
+/*}}}*/
+
+/*Object functions*/
+/*FUNCTION Numericalflux::Configure {{{1*/
+void  Numericalflux::Configure(void* pelementsin,void* pnodesin,void* pmaterialsin){
+
+	DataSet* elementsin=NULL;
+	DataSet* nodesin=NULL;
+
+	/*Recover pointers :*/
+	elementsin=(DataSet*)pelementsin;
+	nodesin=(DataSet*)pnodesin;
+	
+	/*Link this load with its nodes: */
+	if (strcmp(type,"internal")==0){
+		ResolvePointers((Object**)nodes,node_ids,node_offsets,4,nodesin);
+		ResolvePointers((Object**)elements,element_ids,element_offsets,2,elementsin);
+	}
+	else if (strcmp(type,"boundary")==0){
+		ResolvePointers((Object**)nodes,node_ids,node_offsets,2,nodesin);
+		ResolvePointers((Object**)elements,element_ids,element_offsets,1,elementsin);
+	}
+	else ISSMERROR("type not supported yet");
+
+}
+/*}}}*/
+/*FUNCTION Numericalflux::CreateKMatrix {{{1*/
+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);
+	
+	if (strcmp(type,"internal")==0){
+		printf("   element_ids=[%i,%i]\n",element_ids[0],element_ids[1]);
+		printf("   element_offsets=[%i,%i]\n",element_offsets[0],element_offsets[1]);
+		for(i=0;i<2;i++){
+			if(elements[i])elements[i]->Echo();
+		}
+		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("   element_ids=[%i,%i]\n",element_ids[0],element_ids[1]);
+		printf("   element_offsets=[%i,%i]\n",element_offsets[0],element_offsets[1]);
+		for(i=0;i<1;i++){
+			if(elements[i])elements[i]->Echo();
+		}
+		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);
+
+	if (strcmp(type,"internal")==0){
+		printf("   element_ids=[%i,%i]\n",element_ids[0],element_ids[1]);
+		printf("   element_offsets=[%i,%i]\n",element_offsets[0],element_offsets[1]);
+		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("   element_ids=%i\n",element_ids[0]);
+		printf("   element_offsets=%i\n",element_offsets[0]);
+		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::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")){
+		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("type not supported yet");
+
+	/*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::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[0]=-0.5*gauss_coord+0.5;
+	l1l4[1]=+0.5*gauss_coord+0.5;
+
+}
+/*}}}*/
+/*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);
+	}
+	else ISSMERROR("type not supported yet");
+}
+/*}}}*/
+/*FUNCTION Numericalflux::PenaltyCreateKMatrixInternal {{{1*/
+void  Numericalflux::PenaltyCreateKMatrixInternal(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
+
+	ISSMERROR("not supported yet");
+
+}
+/*}}}*/
+/*FUNCTION Numericalflux::PenaltyCreateKMatrixBoundary {{{1*/
+void  Numericalflux::PenaltyCreateKMatrixBoundary(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
+
+	ISSMERROR("type not supported yet");
+
+}
+/*}}}*/
+/*FUNCTION Numericalflux::PenaltyCreatePVector{{{1*/
+void  Numericalflux::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
+
+	if (strcmp(type,"internal")==0){
+
+		PenaltyCreatePVectorInternal(pg,inputs,kmax,analysis_type,sub_analysis_type);
+	}
+	else if (strcmp(type,"boundary")==0){
+
+		PenaltyCreatePVectorBoundary(pg,inputs,kmax,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){
+
+	/*Nothing added to PVector*/
+	return;
+
+}
+/*}}}*/
+/*FUNCTION Numericalflux::PenaltyCreatePVectorBoundary{{{1*/
+void  Numericalflux::PenaltyCreatePVectorBoundary(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
+
+	ISSMERROR("type not supported yet");
+
+}
+/*}}}*/
+/*FUNCTION Numericalflux::UpdateFromInputs {{{1*/
+void  Numericalflux::UpdateFromInputs(void* vinputs){
+
+	/*Do nothing...*/
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Numericalflux.h
===================================================================
--- /issm/trunk/src/c/objects/Numericalflux.h	(revision 3359)
+++ /issm/trunk/src/c/objects/Numericalflux.h	(revision 3359)
@@ -0,0 +1,68 @@
+/*!\file Numericalflux.h
+ * \brief: header file for icefront object
+ */
+
+#ifndef _NUMERICALFLUX_H_
+#define _NUMERICALFLUX_H_
+
+#include "./Load.h"
+#include "./Matpar.h"
+#include "./Element.h"
+#include "./Node.h"
+
+#define NUMERICALFLUXSTRING 20    //max string length
+#define MAX_NUMERICALFLUX_NODES 4 //max number of grids for a certain load
+#define MAX_NUMERICALFLUX_ELEMS 2 //max number of elements for a certain load
+
+class Element;
+class Numericalflux: public Load {
+
+	private: 
+		/*Type of numerical flux*/
+		int  id;
+		char type[NUMERICALFLUXSTRING];
+	
+		/*elements: */
+		Element* elements[MAX_NUMERICALFLUX_ELEMS];
+		int      element_ids[MAX_NUMERICALFLUX_ELEMS];
+		int      element_offsets[MAX_NUMERICALFLUX_ELEMS];
+
+		/*nodes: */
+		Node* nodes[MAX_NUMERICALFLUX_NODES];
+		int   node_ids[MAX_NUMERICALFLUX_NODES];
+		int   node_offsets[MAX_NUMERICALFLUX_NODES];
+
+	public:
+
+		Numericalflux();
+		Numericalflux(int numericalflux_id,char numericalflux_type[NUMERICALFLUXSTRING], int numericalflux_node_ids[MAX_NUMERICALFLUX_NODES],int numericalflux_element_ids[MAX_NUMERICALFLUX_ELEMS]);
+		Object* copy();
+		~Numericalflux();
+
+		void  Echo();
+		void  DeepEcho();
+		void  Marshall(char** pmarshalled_dataset);
+		int   MarshallSize();
+		char* GetName();
+		void  Demarshall(char** pmarshalled_dataset);
+		int   Enum();
+		int   GetId(); 
+		void  GetJacobianDeterminant(double* pJdet,double xyz_list[4][3], double gauss_coord);
+		void  GetNodalFunctions(double* l1l4, double gauss_coord);
+		int   MyRank();
+		void  DistributeNumDofs(int* numdofspernode,int analysis_type,int sub_analysis_type);
+		void  Configure(void* elements,void* nodes,void* materials);
+		void  GetDofList(int* doflist,int* pnumberofdofs);
+		void  UpdateFromInputs(void* inputs);
+		
+		void  CreateKMatrix(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  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);
+};
+
+#endif  /* _NUMERICALFLUX_H_ */
Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 3358)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 3359)
@@ -389,22 +389,23 @@
 
 		CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type);
-
 	}
 	else if (analysis_type==PrognosticAnalysisEnum()){
 
 		CreateKMatrixPrognostic( Kgg,inputs,analysis_type,sub_analysis_type);
-
+	}
+	else if (analysis_type==Prognostic2AnalysisEnum()){
+
+		CreateKMatrixPrognostic2(Kgg,inputs,analysis_type,sub_analysis_type);
 	}
 	else if (analysis_type==BalancedthicknessAnalysisEnum()){
 
 		CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type);
-
 	}
 	else if (analysis_type==BalancedvelocitiesAnalysisEnum()){
 
 		CreateKMatrixBalancedvelocities( Kgg,inputs,analysis_type,sub_analysis_type);
-
 	}
 	else{
+
 		ISSMERROR(exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
 	}
@@ -1604,4 +1605,144 @@
 		}
 #endif
+	} // for (ig=0; ig<num_gauss; ig++)
+
+	/*Add Ke_gg to global matrix Kgg: */
+	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
+
+#ifdef _DEBUGELEMENTS_
+	if(my_rank==RANK && id==ELID){ 
+		printf("      Ke_gg erms:\n");
+		for( i=0; i<numdof; i++){
+			for (j=0;j<numdof;j++){
+				printf("%g ",Ke_gg[i][j]);
+			}
+			printf("\n");
+		}
+		printf("      Ke_gg row_indices:\n");
+		for( i=0; i<numdof; i++){
+			printf("%i ",doflist[i]);
+		}
+
+	}
+#endif
+
+cleanup_and_return: 
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&gauss_weights);
+
+}
+/*}}}*/
+/*FUNCTION CreateKMatrixPrognostic2 {{{1*/
+void  Tria::CreateKMatrixPrognostic2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
+
+	/* local declarations */
+	int             i,j;
+
+	/* node data: */
+	const int    numgrids=3;
+	const int    NDOF1=1;
+	const int    numdof=NDOF1*numgrids;
+	double       xyz_list[numgrids][3];
+	int          doflist[numdof];
+	int          numberofdofspernode;
+
+	/* gaussian points: */
+	int     num_gauss,ig;
+	double* first_gauss_area_coord  =  NULL;
+	double* second_gauss_area_coord =  NULL;
+	double* third_gauss_area_coord  =  NULL;
+	double* gauss_weights           =  NULL;
+	double  gauss_weight;
+	double  gauss_l1l2l3[3];
+
+	/* matrices: */
+	double L[numgrids];
+	double B[2][numgrids];
+	double Bprime[2][numgrids];
+	double DL[2][2]={0.0};
+	double DLprime[2][2]={0.0};
+	double DL_scalar;
+	double Ke_gg[numdof][numdof]={0.0};
+	double Ke_gg1[numdof][numdof]={0.0};
+	double Ke_gg2[numdof][numdof]={0.0};
+	double Jdettria;
+
+	/*input parameters for structural analysis (diagnostic): */
+	double  vx_list[numgrids]={0.0};
+	double  vy_list[numgrids]={0.0};
+	double  vx,vy;
+	double  dt;
+	int     dofs[1]={0};
+	int     found;
+
+	ParameterInputs* inputs=NULL;
+
+	/*recover pointers: */
+	inputs=(ParameterInputs*)vinputs;
+
+	/*recover extra inputs from users, at current convergence iteration: */
+	found=inputs->Recover("vx_average",&vx_list[0],1,dofs,numgrids,(void**)nodes);
+	if(!found)ISSMERROR(" could not find vx_average in inputs!");
+	found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes);
+	if(!found)ISSMERROR(" could not find vy_average in inputs!");
+	found=inputs->Recover("dt",&dt);
+	if(!found)ISSMERROR(" could not find dt in inputs!");
+
+	/* Get node coordinates and dof list: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+	/* Start  looping on the number of gaussian points: */
+	for (ig=0; ig<num_gauss; ig++){
+
+		/*Pick up the gaussian point: */
+		gauss_weight=*(gauss_weights+ig);
+		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
+
+		/*Get L matrix: */
+		GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
+
+		DL_scalar=gauss_weight*Jdettria;
+
+		/*  Do the triple product tL*D*L: */
+		TripleMultiply( &L[0],1,numdof,1,
+					&DL_scalar,1,1,0,
+					&L[0],1,numdof,0,
+					&Ke_gg1[0][0],0);
+
+		/*Get B  and B prime matrix: */
+		/*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
+		GetB_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
+		GetBPrime_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
+
+		//Get vx, vy and their derivatives at gauss point
+		GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
+		GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
+
+		DL_scalar=-dt*gauss_weight*Jdettria;
+
+		DLprime[0][0]=DL_scalar*vx;
+		DLprime[1][1]=DL_scalar*vy;
+
+		//Do the triple product tL*D*L. 
+		TripleMultiply( &B[0][0],2,numdof,1,
+					&DLprime[0][0],2,2,0,
+					&Bprime[0][0],2,numdof,0,
+					&Ke_gg2[0][0],0);
+
+		/* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
+		for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg1[i][j];
+		for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg2[i][j];
+
 	} // for (ig=0; ig<num_gauss; ig++)
 
@@ -1828,5 +1969,4 @@
 		
 			CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
-		
 		}
 		else ISSMERROR(exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
@@ -1839,4 +1979,8 @@
 
 		CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
+	}
+	else if (analysis_type==Prognostic2AnalysisEnum()){
+
+		CreatePVectorPrognostic2( pg,inputs,analysis_type,sub_analysis_type);
 	}
 	else if (analysis_type==BalancedthicknessAnalysisEnum()){
@@ -2321,4 +2465,100 @@
 void  Tria::CreatePVectorPrognostic(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){
 
+
+	/* local declarations */
+	int             i,j;
+
+	/* node data: */
+	const int    numgrids=3;
+	const int    NDOF1=1;
+	const int    numdof=NDOF1*numgrids;
+	double       xyz_list[numgrids][3];
+	int          doflist[numdof];
+	int          numberofdofspernode;
+
+	/* gaussian points: */
+	int     num_gauss,ig;
+	double* first_gauss_area_coord  =  NULL;
+	double* second_gauss_area_coord =  NULL;
+	double* third_gauss_area_coord  =  NULL;
+	double* gauss_weights           =  NULL;
+	double  gauss_weight;
+	double  gauss_l1l2l3[3];
+
+	/* matrix */
+	double pe_g[numgrids]={0.0};
+	double L[numgrids];
+	double Jdettria;
+
+	/*input parameters for structural analysis (diagnostic): */
+	double  accumulation_list[numgrids]={0.0};
+	double  accumulation_g;
+	double  melting_list[numgrids]={0.0};
+	double  melting_g;
+	double  thickness_list[numgrids]={0.0};
+	double  thickness_g;
+	double  dt;
+	int     dofs[1]={0};
+	int     found=0;
+
+	ParameterInputs* inputs=NULL;
+
+	/*recover pointers: */
+	inputs=(ParameterInputs*)vinputs;
+
+	/*recover extra inputs from users, at current convergence iteration: */
+	found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes);
+	if(!found)ISSMERROR(" could not find accumulation in inputs!");
+	found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes);
+	if(!found)ISSMERROR(" could not find melting in inputs!");
+	found=inputs->Recover("thickness",&thickness_list[0],1,dofs,numgrids,(void**)nodes);
+	if(!found)ISSMERROR(" could not find thickness in inputs!");
+	found=inputs->Recover("dt",&dt);
+	if(!found)ISSMERROR(" could not find dt in inputs!");
+
+	/* Get node coordinates and dof list: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+	/* Start  looping on the number of gaussian points: */
+	for (ig=0; ig<num_gauss; ig++){
+		/*Pick up the gaussian point: */
+		gauss_weight=*(gauss_weights+ig);
+		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
+
+		/*Get L matrix: */
+		GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
+
+		/* Get accumulation, melting and thickness at gauss point */
+		GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3);
+		GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3);
+		GetParameterValue(&thickness_g, &thickness_list[0],gauss_l1l2l3);
+
+		/* Add value into pe_g: */
+		for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i];
+
+	} // for (ig=0; ig<num_gauss; ig++)
+
+	/*Add pe_g to global matrix Kgg: */
+	VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
+
+cleanup_and_return: 
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&gauss_weights);
+
+}
+/*}}}*/
+/*FUNCTION CreatePVectorPrognostic2 {{{1*/
+void  Tria::CreatePVectorPrognostic2(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){
 
 	/* local declarations */
Index: /issm/trunk/src/c/objects/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Tria.h	(revision 3358)
+++ /issm/trunk/src/c/objects/Tria.h	(revision 3359)
@@ -125,4 +125,6 @@
 		void  CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
 		void  CreatePVectorPrognostic(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);
+		void  CreateKMatrixPrognostic2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
+		void  CreatePVectorPrognostic2(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);
 		void  CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
 		void  CreatePVectorBalancedthickness(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);
Index: /issm/trunk/src/c/objects/objects.h
===================================================================
--- /issm/trunk/src/c/objects/objects.h	(revision 3358)
+++ /issm/trunk/src/c/objects/objects.h	(revision 3359)
@@ -25,4 +25,5 @@
 #include "./Penpair.h"
 #include "./Pengrid.h"
+#include "./Numericalflux.h"
 #include "./Param.h"
 #include "./Element.h" 
Index: /issm/trunk/src/m/classes/@model/model.m
===================================================================
--- /issm/trunk/src/m/classes/@model/model.m	(revision 3358)
+++ /issm/trunk/src/m/classes/@model/model.m	(revision 3359)
@@ -37,4 +37,5 @@
 	md.nodeconnectivity=NaN;
 	md.elementconnectivity=NaN;
+	md.edges=NaN;
 
 	%Initial 2d mesh 
Index: /issm/trunk/src/m/classes/public/bamg.m
===================================================================
--- /issm/trunk/src/m/classes/public/bamg.m	(revision 3358)
+++ /issm/trunk/src/m/classes/public/bamg.m	(revision 3359)
@@ -280,4 +280,5 @@
 md.y=bamgmesh_out.Vertices(:,2);
 md.elements=bamgmesh_out.Triangles(:,1:3);
+md.edges=bamgmesh_out.ElementEdges;
 md.segments=bamgmesh_out.Segments(:,1:3);
 md.segmentmarkers=bamgmesh_out.Segments(:,4);
Index: /issm/trunk/src/m/classes/public/display/displaymesh.m
===================================================================
--- /issm/trunk/src/m/classes/public/display/displaymesh.m	(revision 3358)
+++ /issm/trunk/src/m/classes/public/display/displaymesh.m	(revision 3359)
@@ -33,10 +33,11 @@
 fielddisplay(md,'y','nodes y coordinate');
 fielddisplay(md,'z','nodes z coordinate');
+fielddisplay(md,'edges','edges of the 2d mesh (node1 node2 element1 element2)');
+
+disp(sprintf('\n      Properties:'));
+fielddisplay(md,'type','mesh type');
 fielddisplay(md,'numlayers','number of extrusion layers');
 fielddisplay(md,'extrusionexponent','exponent for extrusion');
 fielddisplay(md,'dof','maximum number of dofs solved');
-
-disp(sprintf('\n      Properties:'));
-fielddisplay(md,'type','mesh type');
 fielddisplay(md,'bamg','Geometry and 2d mesh properties (if generated by Bamg)');
 fielddisplay(md,'penalties','penalties list');
Index: /issm/trunk/src/m/classes/public/marshall.m
===================================================================
--- /issm/trunk/src/m/classes/public/marshall.m	(revision 3358)
+++ /issm/trunk/src/m/classes/public/marshall.m	(revision 3359)
@@ -82,4 +82,5 @@
 
 WriteData(fid,md.pressureload,'Mat','pressureload');
+WriteData(fid,md.edges,'Mat','edges');
 
 WriteData(fid,md.geothermalflux,'Mat','geothermalflux');
Index: /issm/trunk/src/m/enum/NumericalfluxEnum.m
===================================================================
--- /issm/trunk/src/m/enum/NumericalfluxEnum.m	(revision 3359)
+++ /issm/trunk/src/m/enum/NumericalfluxEnum.m	(revision 3359)
@@ -0,0 +1,9 @@
+function macro=NumericalfluxEnum()
+%NUMERICALFLUXENUM - Enum of Numericalflux
+%
+%   file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum
+%
+%   Usage:
+%      macro=NumericalfluxEnum()
+
+macro=435;
Index: /issm/trunk/src/m/solutions/jpl/prognostic2.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/prognostic2.m	(revision 3358)
+++ /issm/trunk/src/m/solutions/jpl/prognostic2.m	(revision 3359)
@@ -13,5 +13,4 @@
 	displaystring(md.verbose,'%s',['reading prognostic2 model data']);
 	md.analysis_type=Prognostic2AnalysisEnum; models.p=CreateFemModel(md);
-	error('STOP here');
 
 	% figure out number of dof: just for information purposes.
@@ -21,5 +20,6 @@
 	displaystring(md.verbose,'\n%s',['setup inputs...']);
 	inputs=inputlist;
-	inputs=add(inputs,'velocity',models.p.parameters.u_g,'doublevec',3,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'vx',models.p.parameters.vx_g,'doublevec',1,models.p.parameters.numberofnodes);
+	inputs=add(inputs,'vy',models.p.parameters.vy_g,'doublevec',1,models.p.parameters.numberofnodes);
 	inputs=add(inputs,'thickness',models.p.parameters.h_g,'doublevec',1,models.p.parameters.numberofnodes);
 	inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.numberofnodes);
@@ -28,5 +28,5 @@
 
 	displaystring(md.verbose,'\n%s',['call computational core:']);
-	results=prognostic2_core(models,inputs,PrognosticAnalysisEnum(),NoneAnalysisEnum());
+	results=prognostic2_core(models,inputs,Prognostic2AnalysisEnum(),NoneAnalysisEnum());
 
 	displaystring(md.verbose,'\n%s',['load results...']);
Index: /issm/trunk/src/m/solutions/jpl/prognostic2_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/prognostic2_core.m	(revision 3358)
+++ /issm/trunk/src/m/solutions/jpl/prognostic2_core.m	(revision 3359)
@@ -12,10 +12,14 @@
 	displaystring(m.parameters.verbose,'\n%s',['depth averaging velocity...']);
 	%Take only the first two dofs of m.parameters.u_g
-	u_g=get(inputs,'velocity',[1 1 0 0]);
-	u_g=FieldDepthAverage(m.elements,m.nodes,m.loads,m.materials,m.parameters,u_g,'velocity');
-	inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m.parameters.numberofnodes);
+	vx_g=get(inputs,'vx',1);
+	vy_g=get(inputs,'vy',1);
+	vx_g=FieldDepthAverage(m.elements,m.nodes,m.loads,m.materials,m.parameters,vx_g,'vx');
+	vy_g=FieldDepthAverage(m.elements,m.nodes,m.loads,m.materials,m.parameters,vy_g,'vy');
+	inputs=add(inputs,'vx_average',vx_g,'doublevec',1,m.parameters.numberofnodes);
+	inputs=add(inputs,'vy_average',vy_g,'doublevec',1,m.parameters.numberofnodes);
 
 	displaystring(m.parameters.verbose,'\n%s',['call computational core:']);
 	results.h_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type);
+	error('STOP here');
 
 	displaystring(m.parameters.verbose,'\n%s',['extrude computed thickness on all layers:']);