Index: /issm/trunk/src/c/MassFluxx/MassFluxx.cpp
===================================================================
--- /issm/trunk/src/c/MassFluxx/MassFluxx.cpp	(revision 2111)
+++ /issm/trunk/src/c/MassFluxx/MassFluxx.cpp	(revision 2112)
@@ -14,5 +14,5 @@
 
 void MassFluxx(double* pmass_flux, DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials,
-		double* segments,int num_segments,double* vx,double* vy,double* vz){
+		double* segments,int num_segments,double* ug){
 
 	int i,j;
@@ -40,5 +40,5 @@
 			if (element->GetId()==element_id){
 				/*We found the element which owns this segment, use it to compute the mass flux: */
-				mass_flux+=element->MassFlux(segments+5*i+0,vx,vy,vz);
+				mass_flux+=element->MassFlux(segments+5*i+0,ug);
 			}
 		}
@@ -46,6 +46,6 @@
 
 	#ifdef _PARALLEL_
-		MPI_Allreduce ( (void*)&mass_flux,(void*)&all_mass_flux,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
-		all_mass_flux=mass_flux;
+	MPI_Allreduce ( (void*)&mass_flux,(void*)&all_mass_flux,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+	mass_flux=all_mass_flux;
 	#endif
 
Index: /issm/trunk/src/c/MassFluxx/MassFluxx.h
===================================================================
--- /issm/trunk/src/c/MassFluxx/MassFluxx.h	(revision 2111)
+++ /issm/trunk/src/c/MassFluxx/MassFluxx.h	(revision 2112)
@@ -10,5 +10,5 @@
 
 /* local prototypes: */
-void MassFluxx(double* pmass_flux, DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials,double* segments,int num_segments,double* vx,double* vy,double* vz);
+void MassFluxx(double* pmass_flux, DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials,double* segments,int num_segments,double* ug);
 
 
Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp	(revision 2111)
+++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp	(revision 2112)
@@ -89,5 +89,4 @@
 		throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s"," analysis_type: ",iomodel->analysis_type," sub_analysis_type: ",iomodel->sub_analysis_type," not supported yet!"));
 	}
-
 	CreateParametersQmu(pparameters,iomodel,iomodel_handle);
 
Index: /issm/trunk/src/c/ModelProcessorx/Qmu/CreateParametersQmu.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Qmu/CreateParametersQmu.cpp	(revision 2111)
+++ /issm/trunk/src/c/ModelProcessorx/Qmu/CreateParametersQmu.cpp	(revision 2112)
@@ -52,5 +52,5 @@
 		mxArray* pfield2=NULL;
 	#endif
-
+	
 	/*recover parameters : */
 	parameters=*pparameters;
@@ -100,4 +100,5 @@
 		/*Deal with variables for qmu iomodeling: */
 		variabledescriptors=(char**)xmalloc(iomodel->numberofvariables*sizeof(char*));
+
 
 		/*Fetch descriptors: logic varies if we are running parallel or serial. In parallel, qmumarshall 
@@ -229,25 +230,34 @@
 
 				#ifdef _PARALLEL_
-					/*Use the element partitioning vector from the iomodel to down select qmu_mass_flux_segments to only segments that are relevant 
-					 * to this cpu: */
-					my_num_qmu_mass_flux_segments=0;
-					for(j=0;j<num_qmu_mass_flux_segments;j++){
-						if (  (*(qmu_mass_flux_segments+5*j+4))   == iomodel->epart[my_rank])my_qmu_mass_flux_segments++;
-					}
-					if(my_num_qmu_mass_flux_segments){
-						my_qmu_mass_flux_segments=(double*)xcalloc(5*my_num_qmu_mass_flux_segments,sizeof(double));
-						second_count=0;
+
+					/*Only if partitioning exist do we care about the segments: */
+					if(iomodel->epart){
+
+						/*Use the element partitioning vector from the iomodel to down select qmu_mass_flux_segments to only segments that are relevant 
+						 * to this cpu: */
+						my_num_qmu_mass_flux_segments=0;
 						for(j=0;j<num_qmu_mass_flux_segments;j++){
-							if (*(qmu_mass_flux_segments+5*j+4)==iomodel->epart[my_rank]){
-								for(k=0;k<5;k++)*(my_qmu_mass_flux_segments+5*second_count+k)=*(qmu_mass_flux_segments+5*j+k);
-								second_count++;
+							if (  iomodel->epart[(int)(*(qmu_mass_flux_segments+5*j+4))-1]   == my_rank)my_num_qmu_mass_flux_segments++;
+						}
+
+					
+						if(my_num_qmu_mass_flux_segments){
+							my_qmu_mass_flux_segments=(double*)xcalloc(5*my_num_qmu_mass_flux_segments,sizeof(double));
+							second_count=0;
+							for(j=0;j<num_qmu_mass_flux_segments;j++){
+								if (iomodel->epart[(int)*(qmu_mass_flux_segments+5*j+4)-1]==my_rank){
+									for(k=0;k<5;k++)*(my_qmu_mass_flux_segments+5*second_count+k)=*(qmu_mass_flux_segments+5*j+k);
+									second_count++;
+								}
 							}
 						}
+
+						count++;
+						param= new Param(count,"qmu_mass_flux_segments",DOUBLEMAT);
+						param->SetDoubleMat(my_qmu_mass_flux_segments,my_num_qmu_mass_flux_segments,5);
+						parameters->AddObject(param);
+
 					}
-
-					count++;
-					param= new Param(count,"qmu_mass_flux_segments",DOUBLEMAT);
-					param->SetDoubleMat(my_qmu_mass_flux_segments,my_num_qmu_mass_flux_segments,5);
-					parameters->AddObject(param);
+					
 				#else
 
@@ -263,5 +273,4 @@
 			}
 		}
-
 
 		/*Free data: */
@@ -289,5 +298,4 @@
 	}
 
-
 	/*Assign output pointer: */
 	*pparameters=parameters;
Index: /issm/trunk/src/c/Qmux/DakotaResponses.cpp
===================================================================
--- /issm/trunk/src/c/Qmux/DakotaResponses.cpp	(revision 2111)
+++ /issm/trunk/src/c/Qmux/DakotaResponses.cpp	(revision 2112)
@@ -11,8 +11,10 @@
 #include "../DataSet/DataSet.h"    
 #include "../shared/shared.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../MassFluxx/MassFluxx.h"
 
 #undef __FUNCT__ 
 #define __FUNCT__ "DakotaResponses"
-void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Model* model,DataSet* results,int analysis_type,int sub_analysis_type){
+void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Model* model,DataSet* results,DataSet* processed_results,int analysis_type,int sub_analysis_type){
 
 	int i,j;
@@ -25,9 +27,10 @@
 	model->FindParam(&numberofnodes,"numberofnodes");
 
+
 	for(i=0;i<numresponses;i++){
 
 		response_descriptor=responses_descriptors[i];
 
-		//'min_vx' 'max_vx' 'max_abs_vx' 'min_vy' 'max_vy' 'max_abs_vy' 'min_vel' 'max_vel'
+		//'min_vx' 'max_vx' 'max_abs_vx' 'min_vy' 'max_vy' 'max_abs_vy' 'min_vel' 'max_vel, mass_flux'
 
 		if(strcmp(response_descriptor,"min_vel")==0){
@@ -35,5 +38,5 @@
 			double min_vel=0;
 		
-			found=results->FindResult((void*)&vel,"vel");
+			found=processed_results->FindResult((void*)&vel,"vel");
 			if(!found)throw ErrorException(__FUNCT__," could not find vel to compute min_vel");
 
@@ -44,4 +47,8 @@
 
 			if(my_rank==0)responses[i]=min_vel;
+			
+			/*Free ressources:*/
+			xfree((void**)&vel);
+
 			
 		}
@@ -50,5 +57,5 @@
 			double max_vel=0;
 
-			found=results->FindResult((void*)&vel,"vel");
+			found=processed_results->FindResult((void*)&vel,"vel");
 			if(!found)throw ErrorException(__FUNCT__," could not find vel to compute max_vel");
 
@@ -58,4 +65,8 @@
 			}
 			if(my_rank==0)responses[i]=max_vel;
+			
+			/*Free ressources:*/
+			xfree((void**)&vel);
+
 		}
 		else if(strcmp(response_descriptor,"min_vx")==0){
@@ -63,5 +74,5 @@
 			double min_vx=0;
 			
-			found=results->FindResult((void*)&vx,"vx");
+			found=processed_results->FindResult((void*)&vx,"vx");
 			if(!found)throw ErrorException(__FUNCT__," could not find vx to compute min_vx");
 
@@ -71,4 +82,8 @@
 			}
 			if(my_rank==0)responses[i]=min_vx;
+			
+			/*Free ressources:*/
+			xfree((void**)&vx);
+
 		}
 		else if(strcmp(response_descriptor,"max_vx")==0){
@@ -76,5 +91,5 @@
 			double max_vx=0;
 			
-			found=results->FindResult((void*)&vx,"vx");
+			found=processed_results->FindResult((void*)&vx,"vx");
 			if(!found)throw ErrorException(__FUNCT__," could not find vx to compute max_vx");
 
@@ -84,4 +99,8 @@
 			}
 			if(my_rank==0)responses[i]=max_vx;
+			
+			/*Free ressources:*/
+			xfree((void**)&vx);
+
 		}
 		else if(strcmp(response_descriptor,"max_abs_vx")==0){
@@ -89,5 +108,5 @@
 			double max_abs_vx=0;
 			
-			found=results->FindResult((void*)&vx,"vx");
+			found=processed_results->FindResult((void*)&vx,"vx");
 			if(!found)throw ErrorException(__FUNCT__," could not find vx to compute max_abs_vx");
 
@@ -97,4 +116,8 @@
 			}
 			if(my_rank==0)responses[i]=max_abs_vx;
+			
+			/*Free ressources:*/
+			xfree((void**)&vx);
+
 		}
 		else if(strcmp(response_descriptor,"min_vy")==0){
@@ -102,5 +125,5 @@
 			double min_vy=0;
 			
-			found=results->FindResult((void*)&vy,"vy");
+			found=processed_results->FindResult((void*)&vy,"vy");
 			if(!found)throw ErrorException(__FUNCT__," could not find vy to compute min_vy");
 
@@ -110,4 +133,8 @@
 			}
 			if(my_rank==0)responses[i]=min_vy;
+			
+			/*Free ressources:*/
+			xfree((void**)&vy);
+
 		}
 		else if(strcmp(response_descriptor,"max_vy")==0){
@@ -115,5 +142,5 @@
 			double max_vy=0;
 			
-			found=results->FindResult((void*)&vy,"vy");
+			found=processed_results->FindResult((void*)&vy,"vy");
 			if(!found)throw ErrorException(__FUNCT__," could not find vy to compute max_vy");
 
@@ -123,4 +150,8 @@
 			}
 			if(my_rank==0)responses[i]=max_vy;
+			
+			/*Free ressources:*/
+			xfree((void**)&vy);
+
 		}
 		else if(strcmp(response_descriptor,"max_abs_vy")==0){
@@ -128,5 +159,5 @@
 			double max_abs_vy=0;
 			
-			found=results->FindResult((void*)&vy,"vy");
+			found=processed_results->FindResult((void*)&vy,"vy");
 			if(!found)throw ErrorException(__FUNCT__," could not find vy to compute max_abs_vy");
 
@@ -136,4 +167,8 @@
 			}
 			if(my_rank==0)responses[i]=max_abs_vy;
+			
+			/*Free ressources:*/
+			xfree((void**)&vy);
+
 		}
 		else if(strcmp(response_descriptor,"min_vz")==0){
@@ -141,5 +176,5 @@
 			double min_vz=0;
 			
-			found=results->FindResult((void*)&vz,"vz");
+			found=processed_results->FindResult((void*)&vz,"vz");
 			if(!found)throw ErrorException(__FUNCT__," could not find vz to compute min_vz");
 
@@ -149,4 +184,8 @@
 			}
 			if(my_rank==0)responses[i]=min_vz;
+			
+			/*Free ressources:*/
+			xfree((void**)&vz);
+
 		}
 		else if(strcmp(response_descriptor,"max_vz")==0){
@@ -154,5 +193,5 @@
 			double max_vz=0;
 			
-			found=results->FindResult((void*)&vz,"vz");
+			found=processed_results->FindResult((void*)&vz,"vz");
 			if(!found)throw ErrorException(__FUNCT__," could not find vz to compute max_vz");
 
@@ -162,4 +201,8 @@
 			}
 			if(my_rank==0)responses[i]=max_vz;
+			
+			/*Free ressources:*/
+			xfree((void**)&vz);
+
 		}
 		else if(strcmp(response_descriptor,"max_abs_vz")==0){
@@ -167,5 +210,5 @@
 			double max_abs_vz=0;
 			
-			found=results->FindResult((void*)&vz,"vz");
+			found=processed_results->FindResult((void*)&vz,"vz");
 			if(!found)throw ErrorException(__FUNCT__," could not find vz to compute max_abs_vz");
 
@@ -175,4 +218,62 @@
 			}
 			if(my_rank==0)responses[i]=max_abs_vz;
+			
+			/*Free ressources:*/
+			xfree((void**)&vz);
+
+		}
+		else if(strcmp(response_descriptor,"mass_flux")==0){
+
+			int isstokes,ismacayealpattyn,ishutter;
+			Vec       ug=NULL;
+			double*   ug_serial=NULL;
+			double    mass_flux=0;
+			double*   segments=NULL;
+			int       num_segments;
+			DataSet*  elements=NULL;
+			DataSet*  nodes=NULL;
+			DataSet*  materials=NULL;
+			DataSet*  parameters=NULL;
+			FemModel* femmodel=NULL;
+			Param*    param=NULL;
+
+			/*retrieve velocities: */
+			found=results->FindResult(&ug,"u_g");
+			if(!found)throw ErrorException(__FUNCT__," could not find velocity to compute mass_flux");
+			VecToMPISerial(&ug_serial,ug);
+		
+			/*retrieve active fem model: */
+			model->FindParam(&isstokes,"isstokes");
+			model->FindParam(&ismacayealpattyn,"ismacayealpattyn");
+			model->FindParam(&ishutter,"ishutter");
+
+			if(isstokes){
+				femmodel=model->GetFormulation(DiagnosticAnalysisEnum(),StokesAnalysisEnum());
+			}
+			if(ismacayealpattyn){
+				femmodel=model->GetFormulation(DiagnosticAnalysisEnum(),HorizAnalysisEnum());
+			}
+			if(ishutter){
+				femmodel=model->GetFormulation(DiagnosticAnalysisEnum(),HutterAnalysisEnum());
+			}
+
+			/*retrieve qmu_mass_flux_segments: */
+			param=(Param*)femmodel->parameters->FindParamObject("qmu_mass_flux_segments");
+			if(!param)throw ErrorException(__FUNCT__," could not find qmu_mass_flux_segments to compute mass_flux");
+			
+			param->GetParameterValue((void*)&segments);
+			num_segments=param->GetM();
+
+			/*call mass flux module: */
+			MassFluxx(&mass_flux,femmodel->elements,femmodel->nodes,femmodel->loads,femmodel->materials,segments,num_segments,ug_serial);
+			
+			if(my_rank==0)responses[i]=mass_flux;
+			
+			if(my_rank==0)printf("mass_flux %g\n",mass_flux);
+
+			/*Free ressources:*/
+			VecFree(&ug);
+			xfree((void**)&ug_serial);
+			xfree((void**)&segments);
 		}
 		else{
Index: /issm/trunk/src/c/Qmux/Qmux.h
===================================================================
--- /issm/trunk/src/c/Qmux/Qmux.h	(revision 2111)
+++ /issm/trunk/src/c/Qmux/Qmux.h	(revision 2112)
@@ -16,5 +16,5 @@
 void Qmux(Model* model,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
 void SpawnCoreParallel(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, Model* model,ParameterInputs* inputs,int analysis_type,int sub_analysis_type,int counter);
-void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Model* model, DataSet* results,int analysis_type,int sub_analysis_type);
+void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Model* model, DataSet* results,DataSet* processed_results,int analysis_type,int sub_analysis_type);
 #endif
 
Index: /issm/trunk/src/c/Qmux/SpawnCoreParallel.cpp
===================================================================
--- /issm/trunk/src/c/Qmux/SpawnCoreParallel.cpp	(revision 2111)
+++ /issm/trunk/src/c/Qmux/SpawnCoreParallel.cpp	(revision 2112)
@@ -42,4 +42,5 @@
 	/*output from core solutions: */
 	DataSet* results=NULL;
+	DataSet* processed_results=NULL;
 
 	char** responses_descriptors=NULL;
@@ -110,4 +111,28 @@
 	#endif
 
+	/*broadcast response descriptors: */
+	MPI_Bcast(&numresponses,1,MPI_INT,0,MPI_COMM_WORLD); 
+	if(my_rank!=0){
+		responses_descriptors=(char**)xmalloc(numresponses*sizeof(char*));
+	}
+	for(i=0;i<numresponses;i++){
+		if(my_rank==0){
+			string=responses_descriptors[i];
+			string_length=(strlen(string)+1)*sizeof(char);
+		}
+		MPI_Bcast(&string_length,1,MPI_INT,0,MPI_COMM_WORLD); 
+		if(my_rank!=0)string=(char*)xmalloc(string_length);
+		MPI_Bcast(string,string_length,MPI_CHAR,0,MPI_COMM_WORLD); 
+		if(my_rank!=0)responses_descriptors[i]=string;
+	}
+
+	#ifdef _ISSM_DEBUG_
+	for(i=0;i<numresponses;i++){
+		PetscSynchronizedPrintf(MPI_COMM_WORLD,"variable descriptor %i: %s value: %g\n",i,responses_descriptors[i],responses[i]);
+		PetscSynchronizedFlush(MPI_COMM_WORLD);
+	}
+	#endif
+
+
 	_printf_("qmu iteration: %i\n",counter);
 
@@ -150,13 +175,13 @@
 	/*Now process the outputs, before computing the dakota responses: */
 	if(debug)_printf_("process results:\n");
-	ProcessResults(&results,model,analysis_type); 
-	
-
+	ProcessResults(&processed_results,results,model,analysis_type); 
+	
 	/*compute responses on cpu 0: dummy for now! */
 	if(debug)_printf_("compute dakota responses:\n");
-	DakotaResponses(responses,responses_descriptors,numresponses,model,results,analysis_type,sub_analysis_type);
+	DakotaResponses(responses,responses_descriptors,numresponses,model,results,processed_results,analysis_type,sub_analysis_type);
 
 	/*Free ressources:*/
 	delete results;
+	delete processed_results;
 
 	//variables only on cpu != 0
Index: /issm/trunk/src/c/objects/Beam.cpp
===================================================================
--- /issm/trunk/src/c/objects/Beam.cpp	(revision 2111)
+++ /issm/trunk/src/c/objects/Beam.cpp	(revision 2112)
@@ -667,5 +667,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Beam::MassFlux"
-double Beam::MassFlux( double* segment,double* vx,double* vy,double* vz){
-	throw ErrorException(__FUNCT__," not supported yet!");
-}
+double Beam::MassFlux( double* segment,double* ug){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
Index: /issm/trunk/src/c/objects/Beam.h
===================================================================
--- /issm/trunk/src/c/objects/Beam.h	(revision 2111)
+++ /issm/trunk/src/c/objects/Beam.h	(revision 2112)
@@ -87,5 +87,5 @@
 		void  GetParameterValue(double* pvalue, double* value_list,double gauss_coord);
 		void  GetJacobianDeterminant(double* pJdet,double* z_list, double gauss_coord);
-		double MassFlux(double* segment,double* vx,double* vy,double* vz);
+		double MassFlux(double* segment,double* ug);
 
 };
Index: /issm/trunk/src/c/objects/Element.h
===================================================================
--- /issm/trunk/src/c/objects/Element.h	(revision 2111)
+++ /issm/trunk/src/c/objects/Element.h	(revision 2112)
@@ -40,5 +40,5 @@
 		virtual double Misfit(void* inputs,int analysis_type,int sub_analysis_type)=0;
 		virtual void  ComputePressure(Vec p_g)=0;
-		virtual double MassFlux(double* segment,double* vx,double* vy,double* vz)=0;
+		virtual double MassFlux(double* segment,double* ug)=0;
 
 		int           Enum();
Index: /issm/trunk/src/c/objects/Param.cpp
===================================================================
--- /issm/trunk/src/c/objects/Param.cpp	(revision 2111)
+++ /issm/trunk/src/c/objects/Param.cpp	(revision 2112)
@@ -130,5 +130,5 @@
 	
 		case DOUBLEMAT:
-			/*printf("   double matrix. size: %i,%i\n",M,N);
+			printf("   double matrix. size: %i,%i\n",M,N);
 			for(i=0;i<M;i++){
 				for(j=0;j<N;j++){
@@ -136,5 +136,5 @@
 				}
 				printf("\n");
-			}*/
+			}
 			break;
 
Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 2111)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 2112)
@@ -4032,5 +4032,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Penta::MassFlux"
-double Penta::MassFlux( double* segment,double* vx,double* vy,double* vz){
+double Penta::MassFlux( double* segment,double* ug){
 	throw ErrorException(__FUNCT__," not supported yet!");
 }
Index: /issm/trunk/src/c/objects/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Penta.h	(revision 2111)
+++ /issm/trunk/src/c/objects/Penta.h	(revision 2112)
@@ -142,5 +142,5 @@
 		void  CreatePVectorMelting( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type);
 		void  GetPhi(double* phi, double*  epsilon, double viscosity);
-		double MassFlux(double* segment,double* vx,double* vy,double* vz);
+		double MassFlux(double* segment,double* ug);
 
 
Index: /issm/trunk/src/c/objects/Sing.cpp
===================================================================
--- /issm/trunk/src/c/objects/Sing.cpp	(revision 2111)
+++ /issm/trunk/src/c/objects/Sing.cpp	(revision 2112)
@@ -505,5 +505,5 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Sing::MassFlux"
-double Sing::MassFlux( double* segment,double* vx,double* vy,double* vz){
-	throw ErrorException(__FUNCT__," not supported yet!");
-}
+double Sing::MassFlux( double* segment,double* ug){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
Index: /issm/trunk/src/c/objects/Sing.h
===================================================================
--- /issm/trunk/src/c/objects/Sing.h	(revision 2111)
+++ /issm/trunk/src/c/objects/Sing.h	(revision 2112)
@@ -79,5 +79,5 @@
 		void  GradjB(_p_Vec*, void*, int,int);
 		double Misfit(void*,int,int);
-		double MassFlux(double* segment,double* vx,double* vy,double* vz);
+		double MassFlux(double* segment,double* ug);
 
 
Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 2111)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 2112)
@@ -3711,11 +3711,13 @@
 #undef __FUNCT__ 
 #define __FUNCT__ "Tria::MassFlux"
-double Tria::MassFlux( double* segment,double* vx,double* vy,double* vz){
+double Tria::MassFlux( double* segment,double* ug){
 
 	int i;
 
 	const int    numgrids=3;
+	const int    numdofs=2;
+	int          numberofdofspernode;
 	double mass_flux=0;
-	int    doflist[3];
+	int    doflist[numgrids*numdofs];
 	double vx_list[3];
 	double vy_list[3];
@@ -3743,8 +3745,8 @@
 
 	/*recover velocity at three element nodes: */
-	this->GetDofList1(&doflist[0]);
+	this->GetDofList(&doflist[0],&numberofdofspernode);
 	for(i=0;i<3;i++){
-		vx_list[i]=vx[doflist[i]];
-		vy_list[i]=vy[doflist[i]];
+		vx_list[i]=ug[doflist[numberofdofspernode*i+0]];
+		vy_list[i]=ug[doflist[numberofdofspernode*i+1]];
 	}
 
@@ -3775,5 +3777,4 @@
 				  (1.0/3.0*(h1-h2)*(vy1-vy2)+1.0/2.0*h2*(vy1-vy2)+1.0/2.0*(h1-h2)*vy2+h2*vy2)*normal[1]
 				);
-
 	return mass_flux;
 }
Index: /issm/trunk/src/c/objects/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Tria.h	(revision 2111)
+++ /issm/trunk/src/c/objects/Tria.h	(revision 2112)
@@ -122,5 +122,5 @@
 		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);
-		double MassFlux(double* segment,double* vx,double* vy,double* vz);
+		double MassFlux(double* segment,double* ug);
 		double GetArea(void);
 		double GetAreaCoordinate(double x, double y, int which_one);
Index: /issm/trunk/src/c/parallel/ControlTemporaryResults.cpp
===================================================================
--- /issm/trunk/src/c/parallel/ControlTemporaryResults.cpp	(revision 2111)
+++ /issm/trunk/src/c/parallel/ControlTemporaryResults.cpp	(revision 2112)
@@ -19,4 +19,5 @@
 	/*output: */
 	DataSet* temporary_results=NULL;
+	DataSet* results=NULL;
 	Result*  result=NULL;
 	char*    outputfilename=NULL;
@@ -69,7 +70,17 @@
 
 	//process results
-	ProcessResults(&temporary_results,model,ControlAnalysisEnum());
+	ProcessResults(&results,temporary_results,model,ControlAnalysisEnum());
 
 	//Write results on disk
-	OutputResults(temporary_results,outputfilename);
+	OutputResults(results,outputfilename);
+	
+	
+	/*Free ressources:*/
+	delete temporary_results;
+	delete results;
+	xfree((void**)&outputfilename);
+	VecFree(&u_g);
+	xfree((void**)&control_type);
+	xfree((void**)&param_g_copy);
+	xfree((void**)&J_copy);
 }
Index: /issm/trunk/src/c/parallel/ProcessResults.cpp
===================================================================
--- /issm/trunk/src/c/parallel/ProcessResults.cpp	(revision 2111)
+++ /issm/trunk/src/c/parallel/ProcessResults.cpp	(revision 2112)
@@ -21,5 +21,5 @@
 #include "../shared/shared.h"
 
-void ProcessResults(DataSet** presults,Model* model,int analysis_type){
+void ProcessResults(DataSet** pnewresults, DataSet* results,Model* model,int analysis_type){
 
 	int i,n;
@@ -27,7 +27,4 @@
 	Result* newresult=NULL;
 	
-	/*input: */
-	DataSet* results=NULL;
-
 	/*output: */
 	DataSet* newresults=NULL;
@@ -95,7 +92,4 @@
 
 	int numberofnodes;
-
-	/*recover input results: */
-	results=*presults;
 
 	/*Initialize new results: */
@@ -409,9 +403,5 @@
 	}
 
-	/*Delete results: */
-	delete results;
-
-
 	/*Assign output pointers:*/
-	*presults=newresults;
+	*pnewresults=newresults;
 }
Index: /issm/trunk/src/c/parallel/diagnostic.cpp
===================================================================
--- /issm/trunk/src/c/parallel/diagnostic.cpp	(revision 2111)
+++ /issm/trunk/src/c/parallel/diagnostic.cpp	(revision 2112)
@@ -32,4 +32,5 @@
 	/*Results: */
 	DataSet* results=NULL;
+	DataSet* processed_results=NULL;
 	Result* result=NULL;
 	
@@ -123,5 +124,5 @@
 			
 			_printf_("process results:\n");
-			ProcessResults(&results,model,DiagnosticAnalysisEnum());
+			ProcessResults(&processed_results,results,model,DiagnosticAnalysisEnum());
 		}
 		else{
@@ -135,5 +136,5 @@
 
 			_printf_("process results:\n");
-			ProcessResults(&results,model,ControlAnalysisEnum());
+			ProcessResults(&processed_results,results,model,ControlAnalysisEnum());
 		}
 
@@ -168,4 +169,5 @@
 	delete model;
 	delete results;
+	delete processed_results;
 	delete inputs;
 
Index: /issm/trunk/src/c/parallel/parallel.h
===================================================================
--- /issm/trunk/src/c/parallel/parallel.h	(revision 2111)
+++ /issm/trunk/src/c/parallel/parallel.h	(revision 2112)
@@ -52,5 +52,5 @@
 void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,int analysis_type,int sub_analysis_type);
 //int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femmodel,char* filename);
-void ProcessResults(DataSet** presults,Model* model,int analysis_type);
+void ProcessResults(DataSet** pnewresults, DataSet* results,Model* model,int analysis_type);
 
 #endif
Index: /issm/trunk/src/c/parallel/prognostic.cpp
===================================================================
--- /issm/trunk/src/c/parallel/prognostic.cpp	(revision 2111)
+++ /issm/trunk/src/c/parallel/prognostic.cpp	(revision 2112)
@@ -40,4 +40,5 @@
 	/*Results: */
 	DataSet* results=NULL;
+	DataSet* processedresults=NULL;
 	Result*  result=NULL;
 
@@ -123,8 +124,8 @@
 	
 	_printf_("process results:\n");
-	ProcessResults(&results,model,PrognosticAnalysisEnum());
+	ProcessResults(&processedresults,results,model,PrognosticAnalysisEnum());
 	
 	_printf_("write results to disk:\n");
-	OutputResults(results,outputfilename);
+	OutputResults(processedresults,outputfilename);
 
 	_printf_("write lock file:\n");
@@ -138,5 +139,11 @@
 	/*end module: */
 	MODULEEND();
-	
+
+	/*Free ressources:*/
+	delete processedresults;
+	delete results;
+	delete model;
+	delete inputs;
+
 	return 0; //unix success return;
 }
Index: /issm/trunk/src/c/parallel/steadystate.cpp
===================================================================
--- /issm/trunk/src/c/parallel/steadystate.cpp	(revision 2111)
+++ /issm/trunk/src/c/parallel/steadystate.cpp	(revision 2112)
@@ -35,4 +35,5 @@
 	/*Results: */
 	DataSet* results=NULL;
+	DataSet* processed_results=NULL;
 	Result* result=NULL;
 	
@@ -142,5 +143,5 @@
 
 			_printf_("process results:\n");
-			ProcessResults(&results,model,SteadystateAnalysisEnum());
+			ProcessResults(&processed_results,results,model,SteadystateAnalysisEnum());
 		}
 		else{
@@ -154,5 +155,5 @@
 
 			_printf_("process results:\n");
-			ProcessResults(&results,model,ControlAnalysisEnum());
+			ProcessResults(&processed_results,results,model,ControlAnalysisEnum());
 		}
 
@@ -188,4 +189,5 @@
 	delete model;
 	delete results;
+	delete processed_results;
 	delete inputs;
 
Index: /issm/trunk/src/c/parallel/thermal.cpp
===================================================================
--- /issm/trunk/src/c/parallel/thermal.cpp	(revision 2111)
+++ /issm/trunk/src/c/parallel/thermal.cpp	(revision 2112)
@@ -36,4 +36,5 @@
 	/*Results: */
 	DataSet* results=NULL;
+	DataSet* processed_results=NULL;
 	Result*  result=NULL;
 	
@@ -104,8 +105,8 @@
 			
 		_printf_("process results:\n");
-		ProcessResults(&results,model,ThermalAnalysisEnum());
+		ProcessResults(&processed_results,results,model,ThermalAnalysisEnum());
 		
 		_printf_("write results to disk:\n");
-		OutputResults(results,outputfilename);
+		OutputResults(processed_results,outputfilename);
 	}
 	else{
@@ -136,4 +137,5 @@
 	delete model;
 	delete results;
+	delete processed_results;
 	delete inputs;
 	
Index: /issm/trunk/src/c/parallel/transient.cpp
===================================================================
--- /issm/trunk/src/c/parallel/transient.cpp	(revision 2111)
+++ /issm/trunk/src/c/parallel/transient.cpp	(revision 2112)
@@ -32,4 +32,5 @@
 	/*Results: */
 	DataSet* results=NULL;
+	DataSet* processed_results=NULL;
 	Result*  result=NULL;
 	
@@ -130,5 +131,5 @@
 
 		_printf_("process results:\n");
-		ProcessResults(&results,model,TransientAnalysisEnum());
+		ProcessResults(&processed_results,results,model,TransientAnalysisEnum());
 		
 		_printf_("write results to disk:\n");
@@ -151,4 +152,5 @@
 	/*Free ressources:*/
 	delete results;
+	delete processed_results;
 	delete model;
 	delete inputs;
Index: /issm/trunk/src/m/solutions/cielo/SpawnCore.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/SpawnCore.m	(revision 2111)
+++ /issm/trunk/src/m/solutions/cielo/SpawnCore.m	(revision 2112)
@@ -54,5 +54,5 @@
 
 %process results
-results=processresults(models,results);
+processedresults=processresults(models,results);
 
 %now process the results to get response function values
@@ -60,4 +60,4 @@
 for i=1:numel(responsedescriptors),
 	descriptor=responsedescriptors{i};
-	responses(i)=qmuresponse(models,results,descriptor);
+	responses(i)=qmuresponse(models,results,processedresults,descriptor);
 end
Index: /issm/trunk/src/m/solutions/dakota/qmuresponse.m
===================================================================
--- /issm/trunk/src/m/solutions/dakota/qmuresponse.m	(revision 2111)
+++ /issm/trunk/src/m/solutions/dakota/qmuresponse.m	(revision 2112)
@@ -1,17 +1,17 @@
-function response=qmuresponse(models,results,descriptor)
+function response=qmuresponse(models,results,processedresults,descriptor)
 %QMURESPONSE - compute response function from model results.
 
 if strcmpi(descriptor,'max_vel'),
-	response=max(results.vel);
+	response=max(processedresults.vel);
 elseif strcmpi(descriptor,'min_vel'),
-	response=min(results.vel);
+	response=min(processedresults.vel);
 elseif strcmpi(descriptor,'max_vx'),
-	response=max(results.vx);
+	response=max(processedresults.vx);
 elseif strcmpi(descriptor,'min_vx'),
-	response=min(results.vx);
+	response=min(processedresults.vx);
 elseif strcmpi(descriptor,'max_vy'),
-	response=max(results.vy);
+	response=max(processedresults.vy);
 elseif strcmpi(descriptor,'min_vy'),
-	response=min(results.vy);
+	response=min(processedresults.vy);
 elseif strcmpi(descriptor,'mass_flux'),
 	%call mass flux module.
@@ -23,9 +23,9 @@
 	isstokes=m_ds.parameters.isstokes;
 	if ishutter,
-		response=MassFlux(m_dhu.elements,m_dhu.nodes,m_dhu.loads,m_dhu.materials,m_dhu.parameters,results);
+		response=MassFlux(m_dhu.elements,m_dhu.nodes,m_dhu.loads,m_dhu.materials,m_dhu.parameters,results.u_g);
 	elseif ismacayealpattyn,
-		response=MassFlux(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,m_dh.parameters,results);
+		response=MassFlux(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,m_dh.parameters,results.u_g);
 	elseif isstokes,
-		response=MassFlux(m_ds.elements,m_ds.nodes,m_ds.loads,m_ds.materials,m_ds.parameters,results);
+		response=MassFlux(m_ds.elements,m_ds.nodes,m_ds.loads,m_ds.materials,m_ds.parameters,results.u_g);
 	else
 		error('qmuresponse error message: unsupported analysis type for mass_flux computation!');
Index: /issm/trunk/src/mex/MassFlux/MassFlux.cpp
===================================================================
--- /issm/trunk/src/mex/MassFlux/MassFlux.cpp	(revision 2111)
+++ /issm/trunk/src/mex/MassFlux/MassFlux.cpp	(revision 2112)
@@ -17,7 +17,5 @@
 	double*  segments=NULL;
 	int      num_segments;
-	double*  vx=NULL;
-	double*  vy=NULL;
-	double*  vz=NULL;
+	double*  ug=NULL;
 	mxArray* pfield=NULL;
 
@@ -41,11 +39,8 @@
 
 	/*results: */
-	FetchData((void**)&vx,NULL,NULL,mxGetField(RESULTS,0,"vx"),"Matrix","Mat");
-	FetchData((void**)&vy,NULL,NULL,mxGetField(RESULTS,0,"vy"),"Matrix","Mat");
-	pfield=mxGetField(RESULTS,0,"vz");
-	if(pfield)FetchData((void**)&vz,NULL,NULL,pfield,"Matrix","Mat");
+	FetchData((void**)&ug,NULL,NULL,UG,"Matrix","Mat");
 
 	/*!Compute mass flux along the profile: */
-	MassFluxx(&mass_flux, elements,nodes,loads,materials,segments,num_segments,vx,vy,vz);
+	MassFluxx(&mass_flux, elements,nodes,loads,materials,segments,num_segments,ug);
 
 	/*write output datasets: */
@@ -57,7 +52,5 @@
 	delete loads;
 	delete materials;
-	xfree((void**)&vx);
-	xfree((void**)&vy);
-	xfree((void**)&vz);
+	xfree((void**)&ug);
 
 	/*end module: */
Index: /issm/trunk/src/mex/MassFlux/MassFlux.h
===================================================================
--- /issm/trunk/src/mex/MassFlux/MassFlux.h	(revision 2111)
+++ /issm/trunk/src/mex/MassFlux/MassFlux.h	(revision 2112)
@@ -22,5 +22,5 @@
 #define MATERIALS (mxArray*)prhs[3]
 #define PARAMETERS (mxArray*)prhs[4]
-#define RESULTS (mxArray*)prhs[5]
+#define UG (mxArray*)prhs[5]
 
 /* serial output macros: */