Index: /issm/trunk/config.h.in
===================================================================
--- /issm/trunk/config.h.in	(revision 2109)
+++ /issm/trunk/config.h.in	(revision 2110)
@@ -39,7 +39,4 @@
 /* Define to 1 if the X Window System is missing or not being used. */
 #undef X_DISPLAY_MISSING
-
-/* Macro to enable debugging in Dakota. */
-#undef _DEBUG_
 
 /* with Blacs in ISSM src */
@@ -90,2 +87,5 @@
 /* with Triangle in ISSM src */
 #undef _HAVE_TRIANGLE_
+
+/* Macro to enable debugging in ISSM. */
+#undef _ISSM_DEBUG_
Index: /issm/trunk/etc/cluster.rc
===================================================================
--- /issm/trunk/etc/cluster.rc	(revision 2109)
+++ /issm/trunk/etc/cluster.rc	(revision 2110)
@@ -43,5 +43,5 @@
 
 cluster_name=larsen
-cluster_codepath=/u/astrid1/larour/issm/trunk/bin
+cluster_codepath=/u/astrid-r1b/larour/issm/trunk/bin
 cluster_executionpath=/u/wilkes-r1b/larour/Testing/Execution
 
Index: /issm/trunk/scripts/rmlink.sh
===================================================================
--- /issm/trunk/scripts/rmlink.sh	(revision 2110)
+++ /issm/trunk/scripts/rmlink.sh	(revision 2110)
@@ -0,0 +1,8 @@
+#!/bin/bash
+#Erase test decks symbolink links.
+
+for i in `ls | grep _sym | grep -v rmsym | grep -v symlink`
+do
+	rm $i
+done
+
Index: /issm/trunk/scripts/symlink.sh
===================================================================
--- /issm/trunk/scripts/symlink.sh	(revision 2110)
+++ /issm/trunk/scripts/symlink.sh	(revision 2110)
@@ -0,0 +1,8 @@
+#!/bin/bash
+
+#symlinks for test decks, in different categories.
+for i in `ls | grep -v sym | grep -v rmlin`
+do 
+	a=`echo $i | sed 's/_/ /g' | awk '{printf("%s\n",$2);}'`
+	ln -s $i "${a}_sym"
+done
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 2109)
+++ /issm/trunk/src/c/Makefile.am	(revision 2110)
@@ -274,4 +274,6 @@
 					./ElementConnectivityx/ElementConnectivityx.cpp\
 					./ElementConnectivityx/ElementConnectivityx.h\
+					./MassFluxx/MassFluxx.cpp\
+					./MassFluxx/MassFluxx.h\
 					./SystemMatricesx/SystemMatricesx.cpp\
 					./SystemMatricesx/SystemMatricesx.h\
@@ -562,4 +564,6 @@
 					./NormalizeConstraintsx/NormalizeConstraintsx.cpp\
 					./NormalizeConstraintsx/NormalizeConstraintsx.h\
+					./MassFluxx/MassFluxx.cpp\
+					./MassFluxx/MassFluxx.h\
 					./SystemMatricesx/SystemMatricesx.cpp\
 					./SystemMatricesx/SystemMatricesx.h\
Index: /issm/trunk/src/c/MassFluxx/MassFluxx.cpp
===================================================================
--- /issm/trunk/src/c/MassFluxx/MassFluxx.cpp	(revision 2110)
+++ /issm/trunk/src/c/MassFluxx/MassFluxx.cpp	(revision 2110)
@@ -0,0 +1,54 @@
+/*!\file MassFluxx
+ * \brief: compute mass flux along a profile.
+ */
+
+#include "./MassFluxx.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__ "MassFluxx"
+
+#include "../shared/shared.h"
+#include "../include/macros.h"
+#include "../toolkits/toolkits.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+void MassFluxx(double* pmass_flux, DataSet* elements,DataSet* nodes,DataSet* loads,DataSet* materials,
+		double* segments,int num_segments,double* vx,double* vy,double* vz){
+
+	int i,j;
+	extern int num_procs;
+	extern int my_rank;
+	
+	Element* element=NULL;
+	int element_id;
+	
+	/*output: */
+	double mass_flux=0;
+	double all_mass_flux=0;
+
+	/*First, get elements and loads configured: */
+	elements->Configure(elements,loads, nodes, materials);
+	loads->Configure(elements, loads, nodes, materials);
+
+	/*Go through segments, and then elements, and figure out which elements belong to a segment. 
+	 * When we find one, use the element to compute the mass flux on the segment: */
+
+	for(i=0;i<num_segments;i++){
+		element_id=(int)*(segments+5*i+4);
+		for(j=0;j<elements->Size();j++){
+			element=(Element*)elements->GetObjectByOffset(j);
+			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);
+			}
+		}
+	}
+
+	#ifdef _PARALLEL_
+		MPI_Allreduce ( (void*)&mass_flux,(void*)&all_mass_flux,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+		all_mass_flux=mass_flux;
+	#endif
+
+	/*Assign output pointers: */
+	*pmass_flux=mass_flux;
+}
Index: /issm/trunk/src/c/MassFluxx/MassFluxx.h
===================================================================
--- /issm/trunk/src/c/MassFluxx/MassFluxx.h	(revision 2110)
+++ /issm/trunk/src/c/MassFluxx/MassFluxx.h	(revision 2110)
@@ -0,0 +1,16 @@
+/*!\file:  MassFluxx.h
+ * \brief header file for degree of freedoms distribution routines.
+ */ 
+
+#ifndef _MASSFLUXX_H
+#define _MASSFLUXX_H
+
+#include "../DataSet/DataSet.h"
+#include "../objects/objects.h"
+
+/* 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);
+
+
+#endif  /* _MASSFLUXX_H */
+
Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp	(revision 2109)
+++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp	(revision 2110)
@@ -19,8 +19,6 @@
 void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
 
-
 	/*create parameters common to all solutions: */
 	CreateParameters(pparameters,iomodel,iomodel_handle);
-	CreateParametersQmu(pparameters,iomodel,iomodel_handle);
 	CreateParametersControl(pparameters,iomodel,iomodel_handle);
 
@@ -29,10 +27,10 @@
 
 		if (iomodel->sub_analysis_type==HorizAnalysisEnum()){
-			
+
 			CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
 			CreateConstraintsDiagnosticHoriz(pconstraints,iomodel,iomodel_handle);
 			CreateLoadsDiagnosticHoriz(ploads,iomodel,iomodel_handle);
 			CreateParametersDiagnosticHoriz(pparameters,iomodel,iomodel_handle);
-				
+
 		}
 		else if (iomodel->sub_analysis_type==VertAnalysisEnum()){
@@ -91,4 +89,6 @@
 		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 2109)
+++ /issm/trunk/src/c/ModelProcessorx/Qmu/CreateParametersQmu.cpp	(revision 2110)
@@ -17,9 +17,10 @@
 void CreateParametersQmu(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
 	
-	int i;
+	int i,j,k;
 	
 	DataSet* parameters = NULL;
 	Param*   param = NULL;
 	int      count;
+	int      second_count;
 	
 	int*     epart=NULL;
@@ -38,4 +39,12 @@
 	char* qmuerrname=NULL;
 	char* qmuoutname=NULL;
+	extern int my_rank;
+				
+	
+	/*parameters for mass flux: */
+	double* qmu_mass_flux_segments=NULL;
+	double* my_qmu_mass_flux_segments=NULL;
+	int num_qmu_mass_flux_segments=0;
+	int my_num_qmu_mass_flux_segments=0;
 
 	#ifdef _SERIAL_
@@ -57,4 +66,5 @@
 
 	if(iomodel->qmu_analysis){
+
 		//name of qmu input, error and output files
 		qmuinname=(char*)xmalloc((strlen(iomodel->name)+strlen(".qmu.in")+1)*sizeof(char));
@@ -111,4 +121,5 @@
 		#endif
 
+
 		/*Ok, we have all the variable descriptors. Build a parameter with it: */
 		count++;
@@ -140,4 +151,5 @@
 		#endif
 
+
 		/*Ok, we have all the response descriptors. Build a parameter with it: */
 		count++;
@@ -156,4 +168,5 @@
 		#endif
 
+
 		/*partition grids in iomodel->qmu_npart parts, unless a partition is already present: */
 		IoModelFetchData((void**)&dpart,NULL,NULL,iomodel_handle,"part","Matrix","Mat");
@@ -180,4 +193,5 @@
 		param->SetDoubleVec(dpart,iomodel->numberofnodes,1);
 		parameters->AddObject(param);
+
 
 		/*Ok, now if any of the variables input from Dakota are distributed, we are going to need the parameters: */
@@ -204,4 +218,50 @@
 			}
 		}
+
+
+		/*Deal with data needed for some responses: */
+		for(i=0;i<iomodel->numberofresponses;i++){
+			char* descriptor=responsedescriptors[i];
+			if (strcmp(descriptor,"mass_flux")==0){
+
+				/*We need the qmu_mass_flux_segments to be able to compute the mass flux: */
+				IoModelFetchData((void**)&qmu_mass_flux_segments,&num_qmu_mass_flux_segments,NULL,iomodel_handle,"qmu_mass_flux_segments","Matrix","Mat");
+
+				#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;
+						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++;
+							}
+						}
+					}
+
+					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
+
+					count++;
+					param= new Param(count,"qmu_mass_flux_segments",DOUBLEMAT);
+					param->SetDoubleMat(qmu_mass_flux_segments,num_qmu_mass_flux_segments,5);
+					parameters->AddObject(param);
+
+				#endif
+
+				xfree((void**)&qmu_mass_flux_segments);
+				xfree((void**)&my_qmu_mass_flux_segments);
+			}
+		}
+
 
 		/*Free data: */
Index: /issm/trunk/src/c/io/WriteParams.cpp
===================================================================
--- /issm/trunk/src/c/io/WriteParams.cpp	(revision 2109)
+++ /issm/trunk/src/c/io/WriteParams.cpp	(revision 2110)
@@ -28,4 +28,5 @@
 	mwSize		ndim=2;
 	mxArray*    pfield=NULL;
+	mxArray*    pfield2=NULL;
 
 	/*intermediary: */
@@ -109,8 +110,10 @@
 				N=param->GetN();
 				pfield=mxCreateDoubleMatrix(0,0,mxREAL);
-				mxSetM(pfield,M);
-				mxSetN(pfield,N);
+				mxSetM(pfield,N);
+				mxSetN(pfield,M);
 				mxSetPr(pfield,doublemat);
-				mxSetField( dataref, 0, param->GetParameterName(),pfield);
+				//transpose the matrix, written directly to matlab! from C to matlab.
+				mexCallMATLAB(1,&pfield2, 1, &pfield, "'");
+				mxSetField( dataref, 0, param->GetParameterName(),pfield2);
 				break;
 		
Index: /issm/trunk/src/c/issm.h
===================================================================
--- /issm/trunk/src/c/issm.h	(revision 2109)
+++ /issm/trunk/src/c/issm.h	(revision 2110)
@@ -59,4 +59,5 @@
 #include "./ElementConnectivityx/ElementConnectivityx.h"
 #include "./OutputRiftsx/OutputRiftsx.h"
+#include "./MassFluxx/MassFluxx.h"
 
 
Index: /issm/trunk/src/c/objects/Beam.cpp
===================================================================
--- /issm/trunk/src/c/objects/Beam.cpp	(revision 2109)
+++ /issm/trunk/src/c/objects/Beam.cpp	(revision 2110)
@@ -665,3 +665,7 @@
 
 }
-
+#undef __FUNCT__ 
+#define __FUNCT__ "Beam::MassFlux"
+double Beam::MassFlux( double* segment,double* vx,double* vy,double* vz){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
Index: /issm/trunk/src/c/objects/Beam.h
===================================================================
--- /issm/trunk/src/c/objects/Beam.h	(revision 2109)
+++ /issm/trunk/src/c/objects/Beam.h	(revision 2110)
@@ -87,4 +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);
 
 };
Index: /issm/trunk/src/c/objects/Element.h
===================================================================
--- /issm/trunk/src/c/objects/Element.h	(revision 2109)
+++ /issm/trunk/src/c/objects/Element.h	(revision 2110)
@@ -40,4 +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;
 
 		int           Enum();
Index: /issm/trunk/src/c/objects/Param.cpp
===================================================================
--- /issm/trunk/src/c/objects/Param.cpp	(revision 2109)
+++ /issm/trunk/src/c/objects/Param.cpp	(revision 2110)
@@ -272,5 +272,5 @@
 			memcpy(marshalled_dataset,&N,sizeof(N));marshalled_dataset+=sizeof(N);
 			if(M*N){
-				memcpy(marshalled_dataset,doublevec,M*N*sizeof(double));
+				memcpy(marshalled_dataset,doublemat,M*N*sizeof(double));
 				marshalled_dataset+=(M*N*sizeof(double));
 			}
@@ -694,2 +694,17 @@
 
 }
+		
+
+#undef __FUNCT__
+#define __FUNCT__ "SetDoubleMat"
+void  Param::SetDoubleMat(double* value,int pM,int pN){
+	
+	if (type!=DOUBLEMAT) throw ErrorException(__FUNCT__,exprintf("%s%i"," trying to set doublematrix type",type));
+
+	this->M=pM;
+	this->N=pN;
+	if(this->M*this->N){
+		xfree((void**)&doublemat); doublemat=(double*)xcalloc(M*N,sizeof(double));
+		memcpy(doublemat,value,M*N*sizeof(double));
+	}
+}
Index: /issm/trunk/src/c/objects/Param.h
===================================================================
--- /issm/trunk/src/c/objects/Param.h	(revision 2109)
+++ /issm/trunk/src/c/objects/Param.h	(revision 2110)
@@ -50,4 +50,5 @@
 		void  SetDoubleVec(double* value,int size);
 		void  SetDoubleVec(double* value,int size,int ndof);
+		void  SetDoubleMat(double* value,int M,int N);
 		void  SetVec(Vec value);
 		void  SetInteger(int value);
Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 2109)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 2110)
@@ -4028,2 +4028,9 @@
 	*phi=2*pow(epsilon_eff,2.0)*viscosity;
 }
+
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Penta::MassFlux"
+double Penta::MassFlux( double* segment,double* vx,double* vy,double* vz){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
Index: /issm/trunk/src/c/objects/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Penta.h	(revision 2109)
+++ /issm/trunk/src/c/objects/Penta.h	(revision 2110)
@@ -142,4 +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);
 
 
Index: /issm/trunk/src/c/objects/Sing.cpp
===================================================================
--- /issm/trunk/src/c/objects/Sing.cpp	(revision 2109)
+++ /issm/trunk/src/c/objects/Sing.cpp	(revision 2110)
@@ -502,2 +502,8 @@
 
 }
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Sing::MassFlux"
+double Sing::MassFlux( double* segment,double* vx,double* vy,double* vz){
+	throw ErrorException(__FUNCT__," not supported yet!");
+}
Index: /issm/trunk/src/c/objects/Sing.h
===================================================================
--- /issm/trunk/src/c/objects/Sing.h	(revision 2109)
+++ /issm/trunk/src/c/objects/Sing.h	(revision 2110)
@@ -79,4 +79,5 @@
 		void  GradjB(_p_Vec*, void*, int,int);
 		double Misfit(void*,int,int);
+		double MassFlux(double* segment,double* vx,double* vy,double* vz);
 
 
Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 2109)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 2110)
@@ -3709,3 +3709,124 @@
 }
 
-
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::MassFlux"
+double Tria::MassFlux( double* segment,double* vx,double* vy,double* vz){
+
+	int i;
+
+	const int    numgrids=3;
+	double mass_flux=0;
+	int    doflist[3];
+	double vx_list[3];
+	double vy_list[3];
+	double xyz_list[numgrids][3];
+	double gauss_1[3];
+	double gauss_2[3];
+	double normal[2];
+	double length;
+	double x1,y1,x2,y2;
+	double h1,h2;
+	double vx1,vx2,vy1,vy2;
+	double rho_ice;
+	
+	/*Get material parameters :*/
+	rho_ice=this->matpar->GetRhoIce();
+
+	/*First off, check that this segment belongs to this element: */
+	if ((int)*(segment+4)!=this->id)throw ErrorException(__FUNCT__,exprintf("%s%i%s%i","error message: segment with id ",(int)*(segment+4)," does not belong to element with id:",this->id));
+
+	/*Recover segment node locations: */
+	x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
+	
+	/*Get xyz list: */
+	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+
+	/*recover velocity at three element nodes: */
+	this->GetDofList1(&doflist[0]);
+	for(i=0;i<3;i++){
+		vx_list[i]=vx[doflist[i]];
+		vy_list[i]=vy[doflist[i]];
+	}
+
+	/*get area coordinates of 0 and 1 locations: */
+	for(i=0;i<3;i++){
+		gauss_1[i]=this->GetAreaCoordinate(x1,y1,i+1);
+		gauss_2[i]=this->GetAreaCoordinate(x2,y2,i+1);
+	}
+
+	/*get normal of segment: */
+	normal[0]=cos(atan2(x1-x2,y2-y1));
+	normal[1]=sin(atan2(x1-x2,y2-y1));
+
+	/*get length of segment: */
+	length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2));
+
+	/*get thickness and velocity at two segment extremities: */
+	GetParameterValue(&h1, &h[0],gauss_1);
+	GetParameterValue(&h2, &h[0],gauss_2);
+	GetParameterValue(&vx1, &vx_list[0],gauss_1);
+	GetParameterValue(&vy1, &vy_list[0],gauss_1);
+	GetParameterValue(&vx2, &vx_list[0],gauss_2);
+	GetParameterValue(&vy2, &vy_list[0],gauss_2);
+
+
+	mass_flux= rho_ice*length*(  
+				  (1.0/3.0*(h1-h2)*(vx1-vx2)+1.0/2.0*h2*(vx1-vx2)+1.0/2.0*(h1-h2)*vx2+h2*vx2)*normal[0]+
+				  (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;
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::GetArea"
+double Tria::GetArea(void){
+
+	double area=0;
+	const int    numgrids=3;
+	double xyz_list[numgrids][3];
+	double x1,y1,x2,y2,x3,y3;
+
+	/*Get xyz list: */
+	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	x1=xyz_list[0][0]; y1=xyz_list[0][1];
+	x2=xyz_list[1][0]; y2=xyz_list[1][1];
+	x3=xyz_list[2][0]; y3=xyz_list[2][1];
+ 
+	return x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1;
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::GetAreaCoordinate"
+double Tria::GetAreaCoordinate(double x, double y, int which_one){
+
+	double area=0;
+	const int    numgrids=3;
+	double xyz_list[numgrids][3];
+	double x1,y1,x2,y2,x3,y3;
+
+	/*Get area: */
+	area=this->GetArea();
+
+	/*Get xyz list: */
+	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+	x1=xyz_list[0][0]; y1=xyz_list[0][1];
+	x2=xyz_list[1][0]; y2=xyz_list[1][1];
+	x3=xyz_list[2][0]; y3=xyz_list[2][1];
+
+	if(which_one==1){
+		/*Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area*/
+		return ((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area;
+	}
+	else if(which_one==2){
+		/*Get second area coordinate = det(x1-x3  x-x3 ; y1-y3   y-y3)/area*/
+		return ((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area;
+	}
+	else if(which_one==3){
+		/*Get third  area coordinate 1-area1-area2: */
+		return 1-((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area -((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area;
+	}
+	else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n"," error message: area coordinate ",which_one," done not exist!"));
+}
+
+
Index: /issm/trunk/src/c/objects/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Tria.h	(revision 2109)
+++ /issm/trunk/src/c/objects/Tria.h	(revision 2110)
@@ -122,4 +122,7 @@
 		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 GetArea(void);
+		double GetAreaCoordinate(double x, double y, int which_one);
 
 
Index: /issm/trunk/src/c/parallel/diagnostic_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/diagnostic_core.cpp	(revision 2109)
+++ /issm/trunk/src/c/parallel/diagnostic_core.cpp	(revision 2110)
@@ -36,7 +36,9 @@
 	Vec slopey=NULL;
 	Vec riftproperties=NULL;
+	double* u_g_initial=NULL;
 
 	/*flags: */
 	int debug=0;
+	int qmu_analysis=0;
 	int dim=-1;
 	int ishutter=0;
@@ -68,4 +70,5 @@
 	model->FindParam(&stokesreconditioning,"stokesreconditioning");
 	model->FindParam(&numrifts,"numrifts");
+	model->FindParam(&qmu_analysis,"qmu_analysis");
 
 	/*recover fem models: */
@@ -80,4 +83,10 @@
 	fem_sl->FindParam((void*)&numberofdofspernode_sl,"numberofdofspernode");
 	fem_ds->FindParam((void*)&numberofdofspernode_ds,"numberofdofspernode");
+
+	//for qmu analysis, be sure the velocity input we are starting from  is the one in the parameters: */
+	if(qmu_analysis){
+		model->FindParam(&u_g_initial,"u_g",DiagnosticAnalysisEnum(),HorizAnalysisEnum());
+		inputs->Add("velocity",u_g_initial,3,numberofnodes);
+	}
 
 	if(ishutter){
@@ -202,4 +211,5 @@
 	VecFree(&pg);
 	xfree((void**)&dofset);
+	xfree((void**)&u_g_initial);
 
 }
Index: /issm/trunk/src/m/classes/@model/model.m
===================================================================
--- /issm/trunk/src/m/classes/@model/model.m	(revision 2109)
+++ /issm/trunk/src/m/classes/@model/model.m	(revision 2110)
@@ -274,4 +274,6 @@
 	md.variabledescriptors=NaN;
 	md.responsedescriptors=NaN;
+	md.qmu_mass_flux_profile=NaN;
+	md.qmu_mass_flux_segments=NaN;
 
 	%Ice solver string
Index: /issm/trunk/src/m/classes/public/WriteData.m
===================================================================
--- /issm/trunk/src/m/classes/public/WriteData.m	(revision 2109)
+++ /issm/trunk/src/m/classes/public/WriteData.m	(revision 2110)
@@ -33,5 +33,5 @@
 	fwrite(fid,s(2),'int'); 
 	if s(1)*s(2),
-		fwrite(fid,data','double'); %get to the "c" convenction, hence the transpose
+		fwrite(fid,data','double'); %get to the "c" convention, hence the transpose
 	end
 elseif strcmpi(data_type,'Integer'),
Index: /issm/trunk/src/m/classes/public/importancefactors.m
===================================================================
--- /issm/trunk/src/m/classes/public/importancefactors.m	(revision 2109)
+++ /issm/trunk/src/m/classes/public/importancefactors.m	(revision 2110)
@@ -56,2 +56,15 @@
 	factors=importancefactors(npart)';
 end
+
+if numel(factors)==md.numberofgrids,
+	%get areas for each vertex.
+	aire=GetAreas(md.elements,md.x,md.y);
+	num_elements_by_node=md.nodeconnectivity(:,end);
+	grid_aire=zeros(md.numberofgrids,1);
+	for i=1:md.numberofgrids,
+		for j=1:num_elements_by_node(i),
+			grid_aire(i)=grid_aire(i)+aire(md.nodeconnectivity(i,j));
+		end
+	end
+	factors=factors./grid_aire;
+end
Index: /issm/trunk/src/m/classes/public/presolve.m
===================================================================
--- /issm/trunk/src/m/classes/public/presolve.m	(revision 2109)
+++ /issm/trunk/src/m/classes/public/presolve.m	(revision 2110)
@@ -29,2 +29,4 @@
 	count=count+numpairsforthisrift;
 end
+
+
Index: /issm/trunk/src/m/solutions/cielo/CreateFemModel.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/CreateFemModel.m	(revision 2109)
+++ /issm/trunk/src/m/solutions/cielo/CreateFemModel.m	(revision 2110)
@@ -7,8 +7,8 @@
 
 function  m=CreateFEMModel(md)
-
+	
 	displaystring(md.debug,'\n   reading data from model %s...',md.name);
 	[m.elements,m.nodes,m.constraints,m.loads,m.materials,parameters]=ModelProcessor(md);
-	
+
 	displaystring(md.debug,'%s','   generating degrees of freedom...');
 	[m.nodes,m.part,m.tpart]=Dof(m.elements,m.nodes,parameters);
Index: /issm/trunk/src/m/solutions/dakota/preqmu.m
===================================================================
--- /issm/trunk/src/m/solutions/dakota/preqmu.m	(revision 2109)
+++ /issm/trunk/src/m/solutions/dakota/preqmu.m	(revision 2110)
@@ -114,2 +114,5 @@
 md.responsedescriptors=responsedescriptors;
 
+%now, we have to provide all the info necessary for the solutions to compute the responses. For ex, if mass_flux 
+%is a response, we need a profile of points. For max_vel, we don't need anything.
+md=process_qmu_response_data(md);
Index: /issm/trunk/src/m/solutions/dakota/process_qmu_response_data.m
===================================================================
--- /issm/trunk/src/m/solutions/dakota/process_qmu_response_data.m	(revision 2110)
+++ /issm/trunk/src/m/solutions/dakota/process_qmu_response_data.m	(revision 2110)
@@ -0,0 +1,22 @@
+function md=process_qmu_response_data(md)
+%PROCESS_QMU_RESPONSE_DATA - process any data necessary for the solutions to process the data. 
+%
+% Usage: md=process_qmu_response_data(md)
+%
+% See also PREQMU, PRESOLVE
+
+
+
+for i=1:numel(md.responsedescriptors),
+	if strcmpi(md.responsedescriptors{i},'mass_flux'),
+		%we need a profile of points on which to compute the mass_flux, is it here? 
+		if isnans(md.qmu_mass_flux_profile),
+			error('process_qmu_response_data error message: could not find a mass_flux exp profile!');
+		end
+		if ~ischar(md.qmu_mass_flux_profile),
+			error('process_qmu_response_data error message: mass_flux exp profile should be a domain outline name');
+		end
+		%ok, process the qmu_mass_flux_profile to build a list of segments: 
+		md.qmu_mass_flux_segments=MassFluxProcessProfile(md);
+	end
+end
Index: /issm/trunk/src/m/solutions/dakota/qmumarshall.m
===================================================================
--- /issm/trunk/src/m/solutions/dakota/qmumarshall.m	(revision 2109)
+++ /issm/trunk/src/m/solutions/dakota/qmumarshall.m	(revision 2110)
@@ -69,4 +69,17 @@
 end
 
+%write response specific data
+count=0;
+for i=1:length(response_fieldnames),
+	field_name=response_fieldnames{i};
+	fieldresponses=responses.(field_name);
+	for j=1:length(fieldresponses),
+		descriptor=fieldresponses(j).descriptor;
+		if strcmpi(descriptor,'mass_flux'),
+			WriteData(fid,md.qmu_mass_flux_segments,'Mat','qmu_mass_flux_segments');
+		end
+	end
+end
+
 %write part and npart to disk
 WriteData(fid,md.npart,'Integer','npart');
Index: /issm/trunk/src/m/solutions/dakota/qmuresponse.m
===================================================================
--- /issm/trunk/src/m/solutions/dakota/qmuresponse.m	(revision 2109)
+++ /issm/trunk/src/m/solutions/dakota/qmuresponse.m	(revision 2110)
@@ -4,4 +4,31 @@
 if strcmpi(descriptor,'max_vel'),
 	response=max(results.vel);
+elseif strcmpi(descriptor,'min_vel'),
+	response=min(results.vel);
+elseif strcmpi(descriptor,'max_vx'),
+	response=max(results.vx);
+elseif strcmpi(descriptor,'min_vx'),
+	response=min(results.vx);
+elseif strcmpi(descriptor,'max_vy'),
+	response=max(results.vy);
+elseif strcmpi(descriptor,'min_vy'),
+	response=min(results.vy);
+elseif strcmpi(descriptor,'mass_flux'),
+	%call mass flux module.
+	m_dh=models.dh;
+	m_dhu=models.dhu;
+	m_ds=models.ds;
+	ishutter=m_dh.parameters.ishutter;
+	ismacayealpattyn=m_dhu.parameters.ismacayealpattyn;
+	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);
+	elseif ismacayealpattyn,
+		response=MassFlux(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,m_dh.parameters,results);
+	elseif isstokes,
+		response=MassFlux(m_ds.elements,m_ds.nodes,m_ds.loads,m_ds.materials,m_ds.parameters,results);
+	else
+		error('qmuresponse error message: unsupported analysis type for mass_flux computation!');
+	end
 else
 	error(['qmuresponse error message: unknow descriptor ' descriptor]);
Index: /issm/trunk/src/m/solutions/dakota/responsedataprocessing/MassFluxProcessProfile.m
===================================================================
--- /issm/trunk/src/m/solutions/dakota/responsedataprocessing/MassFluxProcessProfile.m	(revision 2110)
+++ /issm/trunk/src/m/solutions/dakota/responsedataprocessing/MassFluxProcessProfile.m	(revision 2110)
@@ -0,0 +1,14 @@
+function segments=MassFluxProcessProfile(md);
+%MASSFLUXPROCESSPROFILE: process an argus domain outlien profile into a list of segments.
+%
+% Usage: segments=MassFluxProcessProfile(md);
+%
+%
+% See also: PROCESS_QMU_RESPONSE_DATA, PREQMU
+
+
+%first read the profile points.
+profile=expread(['../' md.qmu_mass_flux_profile],1); %go to ../ because the exp file is in the root directory.
+
+%project this profile onto mesh.
+segments=ProfileProjectOntoMesh(md,profile);
Index: /issm/trunk/src/m/utils/Mesh/NodeInElement.m
===================================================================
--- /issm/trunk/src/m/utils/Mesh/NodeInElement.m	(revision 2110)
+++ /issm/trunk/src/m/utils/Mesh/NodeInElement.m	(revision 2110)
@@ -0,0 +1,41 @@
+function node_in_element=NodeInElement(newx,newy,elements,x,y,nodeconnectivity);
+%NODEINELEMENT: find for a list of nodes (in newx,newy), which elements in the mesh (elements,x,y) they belong to.
+%
+%  Usage: node_in_element=NodeInElement(newx,newy,elements,x,y,md.nodeconnectivity);
+%
+%  See also Nodeconnectivity
+%
+epsilon=10^-10;
+
+%compute some quantities that will speed up the process
+x3x1=x(elements(:,1))-x(elements(:,3));
+y3y1=y(elements(:,1))-y(elements(:,3));
+x3x2=x(elements(:,2))-x(elements(:,3));
+y3y2=y(elements(:,2))-y(elements(:,3));
+x3=x(elements(:,3));
+y3=y(elements(:,3));
+delta=x(elements(:,2)).*y(elements(:,3))-y(elements(:,2)).*x(elements(:,3))-x(elements(:,1)).*y(elements(:,3))+y(elements(:,1)).*x(elements(:,3))+x(elements(:,1)).*y(elements(:,2))-y(elements(:,1)).*x(elements(:,2));
+
+%max connectivity:
+max_connectivity=max(nodeconnectivity(:,end));
+node_in_element=zeros(length(newx),max_connectivity+1); %last column is the number of elements to which the row node is connected.
+
+for i=1:length(newx),
+	x0=newx(i);
+	y0=newy(i);
+	
+	%first area coordinate
+	area_1=(y3y2.*(x0-x3)-x3x2.*(y0-y3))./delta;
+	%second area coordinate
+	area_2=(x3x1.*(y0-y3)-y3y1.*(x0-x3))./delta;
+	%third area coordinate
+	area_3=1-area_1-area_2;
+	
+	%get elements for which all area coordinates are positive (meaning (x0,y0) belongs to these elements
+	pos=find((area_1>=0-epsilon) & (area_2>=0-epsilon) & (area_3>=0-epsilon));
+
+	num_elements=length(pos);
+
+	node_in_element(i,1:num_elements)=pos;
+	node_in_element(i,end)=num_elements;
+end
Index: /issm/trunk/src/m/utils/Mesh/ProfileProjectOntoMesh.m
===================================================================
--- /issm/trunk/src/m/utils/Mesh/ProfileProjectOntoMesh.m	(revision 2110)
+++ /issm/trunk/src/m/utils/Mesh/ProfileProjectOntoMesh.m	(revision 2110)
@@ -0,0 +1,52 @@
+function mesh_profile=ProfileProjectOntoMesh(md,profile)
+%PROFILEPROJECTONTOMESH: project a profile (made of arbitrary points) onto a mesh, so that we end 
+%                        up with a list of segments self contained onto elements.
+%
+% Usage: mesh_profile=ProfileProjectOntoMesh(md,profile)
+%
+% See also intersections.m
+
+%make a curve out of the mesh, to use the intersections routine.
+rows=[md.elements md.elements(:,1)]'; rows=rows(:);
+x=md.x(rows);
+y=md.y(rows);
+
+%[x0,y0] = intersections(profile.x,profile.y,x,y,1);
+[x0,y0,indices,j] = intersections(profile.x,profile.y,x,y);
+
+%process x0,y0 so they do not include profile.x or profile.y
+processed_indices=[];
+processed_x=[];
+processed_y=[];
+for i=1:numel(indices),
+	if(((indices(i)-floor(indices(i)))~=0) & ((ceil(indices(i))-indices(i))~=0))
+		processed_indices=[processed_indices;floor(indices(i))];
+		processed_x=[processed_x;x0(i)];
+		processed_y=[processed_y;y0(i)];
+	end
+end
+
+%now merge profile.x,profile.y with processed_x,processed_y, at locatoins processed_indices:
+newx=profile.x;
+newy=profile.y;
+
+count=1;
+for i=1:numel(profile.x),
+	pos=find(processed_indices==i);
+	if ~isempty(pos),
+		newx=[newx(1:count); processed_x(pos); newx(count+1:end)];
+		newy=[newy(1:count); processed_y(pos); newy(count+1:end)];
+		count=count+length(pos)+1;
+	end
+end
+
+%now, for each node, figure out which element it belongs to.
+node_in_element=NodeInElement(newx,newy,md.elements,md.x,md.y,md.nodeconnectivity);
+
+mesh_profile=[newx(1:end-1) newy(1:end-1) newx(2:end) newy(2:end) zeros(length(newy(2:end)),1)];
+
+%find which element each segment belongs to.
+for i=1:length(newx)-1,
+	common=intersect(node_in_element(i,1:node_in_element(i,end)), node_in_element(i+1,1:node_in_element(i+1,end)));
+	mesh_profile(i,end)=common(1);
+end
Index: /issm/trunk/src/m/utils/Mesh/intersections.m
===================================================================
--- /issm/trunk/src/m/utils/Mesh/intersections.m	(revision 2110)
+++ /issm/trunk/src/m/utils/Mesh/intersections.m	(revision 2110)
@@ -0,0 +1,268 @@
+function [x0,y0,iout,jout] = intersections(x1,y1,x2,y2,robust)
+%INTERSECTIONS Intersections of curves.
+%   Computes the (x,y) locations where two curves intersect.  The curves
+%   can be broken with NaNs or have vertical segments.
+%
+% Example:
+%   [X0,Y0] = intersections(X1,Y1,X2,Y2,ROBUST);
+%
+% where X1 and Y1 are equal-length vectors of at least two points and
+% represent curve 1.  Similarly, X2 and Y2 represent curve 2.
+% X0 and Y0 are column vectors containing the points at which the two
+% curves intersect.
+%
+% ROBUST (optional) set to 1 or true means to use a slight variation of the
+% algorithm that might return duplicates of some intersection points, and
+% then remove those duplicates.  The default is true, but since the
+% algorithm is slightly slower you can set it to false if you know that
+% your curves don't intersect at any segment boundaries.  Also, the robust
+% version properly handles parallel and overlapping segments.
+%
+% The algorithm can return two additional vectors that indicate which
+% segment pairs contain intersections and where they are:
+%
+%   [X0,Y0,I,J] = intersections(X1,Y1,X2,Y2,ROBUST);
+%
+% For each element of the vector I, I(k) = (segment number of (X1,Y1)) +
+% (how far along this segment the intersection is).  For example, if I(k) =
+% 45.25 then the intersection lies a quarter of the way between the line
+% segment connecting (X1(45),Y1(45)) and (X1(46),Y1(46)).  Similarly for
+% the vector J and the segments in (X2,Y2).
+%
+% You can also get intersections of a curve with itself.  Simply pass in
+% only one curve, i.e.,
+%
+%   [X0,Y0] = intersections(X1,Y1,ROBUST);
+%
+% where, as before, ROBUST is optional.
+
+% Version: 1.10, 25 February 2008
+% Author:  Douglas M. Schwarz
+% Email:   dmschwarz=ieee*org, dmschwarz=urgrad*rochester*edu
+% Real_email = regexprep(Email,{'=','*'},{'@','.'})
+
+
+% Theory of operation:
+%
+% Given two line segments, L1 and L2,
+%
+%   L1 endpoints:  (x1(1),y1(1)) and (x1(2),y1(2))
+%   L2 endpoints:  (x2(1),y2(1)) and (x2(2),y2(2))
+%
+% we can write four equations with four unknowns and then solve them.  The
+% four unknowns are t1, t2, x0 and y0, where (x0,y0) is the intersection of
+% L1 and L2, t1 is the distance from the starting point of L1 to the
+% intersection relative to the length of L1 and t2 is the distance from the
+% starting point of L2 to the intersection relative to the length of L2.
+%
+% So, the four equations are
+%
+%    (x1(2) - x1(1))*t1 = x0 - x1(1)
+%    (x2(2) - x2(1))*t2 = x0 - x2(1)
+%    (y1(2) - y1(1))*t1 = y0 - y1(1)
+%    (y2(2) - y2(1))*t2 = y0 - y2(1)
+%
+% Rearranging and writing in matrix form,
+%
+%  [x1(2)-x1(1)       0       -1   0;      [t1;      [-x1(1);
+%        0       x2(2)-x2(1)  -1   0;   *   t2;   =   -x2(1);
+%   y1(2)-y1(1)       0        0  -1;       x0;       -y1(1);
+%        0       y2(2)-y2(1)   0  -1]       y0]       -y2(1)]
+%
+% Let's call that A*T = B.  We can solve for T with T = A\B.
+%
+% Once we have our solution we just have to look at t1 and t2 to determine
+% whether L1 and L2 intersect.  If 0 <= t1 < 1 and 0 <= t2 < 1 then the two
+% line segments cross and we can include (x0,y0) in the output.
+%
+% In principle, we have to perform this computation on every pair of line
+% segments in the input data.  This can be quite a large number of pairs so
+% we will reduce it by doing a simple preliminary check to eliminate line
+% segment pairs that could not possibly cross.  The check is to look at the
+% smallest enclosing rectangles (with sides parallel to the axes) for each
+% line segment pair and see if they overlap.  If they do then we have to
+% compute t1 and t2 (via the A\B computation) to see if the line segments
+% cross, but if they don't then the line segments cannot cross.  In a
+% typical application, this technique will eliminate most of the potential
+% line segment pairs.
+
+
+% Input checks.
+error(nargchk(2,5,nargin))
+
+% Adjustments when fewer than five arguments are supplied.
+switch nargin
+	case 2
+		robust = true;
+		x2 = x1;
+		y2 = y1;
+		self_intersect = true;
+	case 3
+		robust = x2;
+		x2 = x1;
+		y2 = y1;
+		self_intersect = true;
+	case 4
+		robust = true;
+		self_intersect = false;
+	case 5
+		self_intersect = false;
+end
+
+% x1 and y1 must be vectors with same number of points (at least 2).
+if sum(size(x1) > 1) ~= 1 || sum(size(y1) > 1) ~= 1 || ...
+		length(x1) ~= length(y1)
+	error('X1 and Y1 must be equal-length vectors of at least 2 points.')
+end
+% x2 and y2 must be vectors with same number of points (at least 2).
+if sum(size(x2) > 1) ~= 1 || sum(size(y2) > 1) ~= 1 || ...
+		length(x2) ~= length(y2)
+	error('X2 and Y2 must be equal-length vectors of at least 2 points.')
+end
+
+
+% Force all inputs to be column vectors.
+x1 = x1(:);
+y1 = y1(:);
+x2 = x2(:);
+y2 = y2(:);
+
+% Compute number of line segments in each curve and some differences we'll
+% need later.
+n1 = length(x1) - 1;
+n2 = length(x2) - 1;
+xy1 = [x1 y1];
+xy2 = [x2 y2];
+dxy1 = diff(xy1);
+dxy2 = diff(xy2);
+
+% Determine the combinations of i and j where the rectangle enclosing the
+% i'th line segment of curve 1 overlaps with the rectangle enclosing the
+% j'th line segment of curve 2.
+[i,j] = find(repmat(min(x1(1:end-1),x1(2:end)),1,n2) <= ...
+	repmat(max(x2(1:end-1),x2(2:end)).',n1,1) & ...
+	repmat(max(x1(1:end-1),x1(2:end)),1,n2) >= ...
+	repmat(min(x2(1:end-1),x2(2:end)).',n1,1) & ...
+	repmat(min(y1(1:end-1),y1(2:end)),1,n2) <= ...
+	repmat(max(y2(1:end-1),y2(2:end)).',n1,1) & ...
+	repmat(max(y1(1:end-1),y1(2:end)),1,n2) >= ...
+	repmat(min(y2(1:end-1),y2(2:end)).',n1,1));
+
+% Find segments pairs which have at least one vertex = NaN and remove them.
+% This line is a fast way of finding such segment pairs.  We take
+% advantage of the fact that NaNs propagate through calculations, in
+% particular subtraction (in the calculation of dxy1 and dxy2, which we
+% need anyway) and addition.
+% At the same time we can remove redundant combinations of i and j in the
+% case of finding intersections of a line with itself.
+if self_intersect
+	remove = isnan(sum(dxy1(i,:) + dxy2(j,:),2)) | j <= i + 1;
+else
+	remove = isnan(sum(dxy1(i,:) + dxy2(j,:),2));
+end
+i(remove) = [];
+j(remove) = [];
+
+% Initialize matrices.  We'll put the T's and B's in matrices and use them
+% one column at a time.  AA is a 3-D extension of A where we'll use one
+% plane at a time.
+n = length(i);
+T = zeros(4,n);
+AA = zeros(4,4,n);
+AA([1 2],3,:) = -1;
+AA([3 4],4,:) = -1;
+AA([1 3],1,:) = dxy1(i,:).';
+AA([2 4],2,:) = dxy2(j,:).';
+B = -[x1(i) x2(j) y1(i) y2(j)].';
+
+% Loop through possibilities.  Trap singularity warning and then use
+% lastwarn to see if that plane of AA is near singular.  Process any such
+% segment pairs to determine if they are colinear (overlap) or merely
+% parallel.  That test consists of checking to see if one of the endpoints
+% of the curve 2 segment lies on the curve 1 segment.  This is done by
+% checking the cross product
+%
+%   (x1(2),y1(2)) - (x1(1),y1(1)) x (x2(2),y2(2)) - (x1(1),y1(1)).
+%
+% If this is close to zero then the segments overlap.
+
+% If the robust option is false then we assume no two segment pairs are
+% parallel and just go ahead and do the computation.  If A is ever singular
+% a warning will appear.  This is faster and obviously you should use it
+% only when you know you will never have overlapping or parallel segment
+% pairs.
+
+if robust
+	overlap = false(1,n);
+	warning_state = warning('off','MATLAB:singularMatrix');
+	% Use try-catch to guarantee original warning state is restored.
+	try
+		lastwarn('')
+		for k = 1:n
+			T(:,k) = AA(:,:,k)\B(:,k);
+			[unused,last_warn] = lastwarn;
+			lastwarn('')
+			if strcmp(last_warn,'MATLAB:singularMatrix')
+				% Force in_range(k) to be false.
+				T(1,k) = NaN;
+				% Determine if these segments overlap or are just parallel.
+				overlap(k) = rcond([dxy1(i(k),:);xy2(j(k),:) - xy1(i(k),:)]) < eps;
+			end
+		end
+		warning(warning_state)
+	catch
+		warning(warning_state)
+		rethrow(lasterror)
+	end
+	% Find where t1 and t2 are between 0 and 1 and return the corresponding
+	% x0 and y0 values.
+	in_range = T(1,:) >= 0 & T(2,:) >= 0 & T(1,:) <= 1 & T(2,:) <= 1;
+	% For overlapping segment pairs the algorithm will return an
+	% intersection point that is at the center of the overlapping region.
+	if any(overlap)
+		ia = i(overlap);
+		ja = j(overlap);
+		% set x0 and y0 to middle of overlapping region.
+		T(3,overlap) = (max(min(x1(ia),x1(ia+1)),min(x2(ja),x2(ja+1))) + ...
+			min(max(x1(ia),x1(ia+1)),max(x2(ja),x2(ja+1)))).'/2;
+		T(4,overlap) = (max(min(y1(ia),y1(ia+1)),min(y2(ja),y2(ja+1))) + ...
+			min(max(y1(ia),y1(ia+1)),max(y2(ja),y2(ja+1)))).'/2;
+		selected = in_range | overlap;
+	else
+		selected = in_range;
+	end
+	xy0 = T(3:4,selected).';
+	
+	% Remove duplicate intersection points.
+	[xy0,index] = unique(xy0,'rows');
+	x0 = xy0(:,1);
+	y0 = xy0(:,2);
+	
+	% Compute how far along each line segment the intersections are.
+	if nargout > 2
+		sel_index = find(selected);
+		sel = sel_index(index);
+		iout = i(sel) + T(1,sel).';
+		jout = j(sel) + T(2,sel).';
+	end
+else % non-robust option
+	for k = 1:n
+		[L,U] = lu(AA(:,:,k));
+		T(:,k) = U\(L\B(:,k));
+	end
+	
+	% Find where t1 and t2 are between 0 and 1 and return the corresponding
+	% x0 and y0 values.
+	in_range = T(1,:) >= 0 & T(2,:) >= 0 & T(1,:) < 1 & T(2,:) < 1;
+	x0 = T(3,in_range).';
+	y0 = T(4,in_range).';
+	
+	% Compute how far along each line segment the intersections are.
+	if nargout > 2
+		iout = i(in_range) + T(1,in_range).';
+		jout = j(in_range) + T(2,in_range).';
+	end
+end
+
+% Plot the results (useful for debugging).
+% plot(x1,y1,x2,y2,x0,y0,'ok');
Index: /issm/trunk/src/mex/Makefile.am
===================================================================
--- /issm/trunk/src/mex/Makefile.am	(revision 2109)
+++ /issm/trunk/src/mex/Makefile.am	(revision 2110)
@@ -22,4 +22,5 @@
 				InterpFromMesh2d \
 				InterpFromMesh3d \
+				MassFlux\
 				Mergesolutionfromftog\
 				MeshPartition\
@@ -137,4 +138,7 @@
 			  Mergesolutionfromftog/Mergesolutionfromftog.h
 
+MassFlux_SOURCES = MassFlux/MassFlux.cpp\
+			  MassFlux/MassFlux.h
+
 MeshPartition_SOURCES = MeshPartition/MeshPartition.cpp\
 			  MeshPartition/MeshPartition.h
Index: /issm/trunk/src/mex/MassFlux/MassFlux.cpp
===================================================================
--- /issm/trunk/src/mex/MassFlux/MassFlux.cpp	(revision 2110)
+++ /issm/trunk/src/mex/MassFlux/MassFlux.cpp	(revision 2110)
@@ -0,0 +1,72 @@
+/*\file MassFlux.c
+ *\brief: compute mass flux along a profile.
+ */
+
+#include "./MassFlux.h"
+
+void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){
+
+	/*diverse: */
+	int   noerr=1;
+
+	/*input datasets: */
+	DataSet* elements=NULL;
+	DataSet* nodes=NULL;
+	DataSet* loads=NULL;
+	DataSet* materials=NULL;
+	double*  segments=NULL;
+	int      num_segments;
+	double*  vx=NULL;
+	double*  vy=NULL;
+	double*  vz=NULL;
+	mxArray* pfield=NULL;
+
+	/* output datasets: */
+	double   mass_flux;
+
+	/*Boot module: */
+	MODULEBOOT();
+
+	/*checks on arguments on the matlab side: */
+	CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&MassFluxUsage);
+
+	/*Input datasets: */
+	FetchData((void**)&elements,NULL,NULL,ELEMENTS,"DataSet",NULL);
+	FetchData((void**)&nodes,NULL,NULL,NODES,"DataSet",NULL);
+	FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL);
+	FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL);
+
+	/*parameters: */
+	FetchData((void**)&segments,&num_segments,NULL,mxGetField(PARAMETERS,0,"qmu_mass_flux_segments"),"Matrix","Mat");
+
+	/*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");
+
+	/*!Compute mass flux along the profile: */
+	MassFluxx(&mass_flux, elements,nodes,loads,materials,segments,num_segments,vx,vy,vz);
+
+	/*write output datasets: */
+	WriteData(RESPONSE,&mass_flux,0,0,"Scalar",NULL);
+	
+	/*Free ressources: */
+	delete elements;
+	delete nodes;
+	delete loads;
+	delete materials;
+	xfree((void**)&vx);
+	xfree((void**)&vy);
+	xfree((void**)&vz);
+
+	/*end module: */
+	MODULEEND();
+}
+
+void MassFluxUsage(void)
+{
+	_printf_("\n");
+	_printf_("   usage: [Kgg,pg] = %s(eleemnts,nodes,loads,materials,params,inputs,analysis_type);\n",__FUNCT__);
+	_printf_("\n");
+}
Index: /issm/trunk/src/mex/MassFlux/MassFlux.h
===================================================================
--- /issm/trunk/src/mex/MassFlux/MassFlux.h	(revision 2110)
+++ /issm/trunk/src/mex/MassFlux/MassFlux.h	(revision 2110)
@@ -0,0 +1,36 @@
+
+/*
+	MassFlux.h
+*/
+
+
+#ifndef _MASSFLUX_H
+#define _MASSFLUX_H
+
+/* local prototypes: */
+void MassFluxUsage(void);
+
+#include "../../c/issm.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__  "MassFlux"
+
+/* serial input macros: */
+#define ELEMENTS (mxArray*)prhs[0]
+#define NODES (mxArray*)prhs[1]
+#define LOADS (mxArray*)prhs[2]
+#define MATERIALS (mxArray*)prhs[3]
+#define PARAMETERS (mxArray*)prhs[4]
+#define RESULTS (mxArray*)prhs[5]
+
+/* serial output macros: */
+#define RESPONSE (mxArray**)&plhs[0]
+
+/* serial arg counts: */
+#undef NLHS
+#define NLHS  1
+#undef NRHS
+#define NRHS  6
+
+
+#endif  /* _MASSFLUX_H */
Index: /issm/trunk/test/Verification/IceSheetIceFrontM2dDakota_25/Square.par
===================================================================
--- /issm/trunk/test/Verification/IceSheetIceFrontM2dDakota_25/Square.par	(revision 2109)
+++ /issm/trunk/test/Verification/IceSheetIceFrontM2dDakota_25/Square.par	(revision 2110)
@@ -60,3 +60,3 @@
 md.qmu_analysis=1;
 
-md.eps_rel=10^-15; %tighten for qmu analysese
+md.eps_rel=10^-10; %tighten for qmu analysese
Index: /issm/trunk/test/Verification/IceSheetIceFrontM2dDakota_25/runme.m
===================================================================
--- /issm/trunk/test/Verification/IceSheetIceFrontM2dDakota_25/runme.m	(revision 2109)
+++ /issm/trunk/test/Verification/IceSheetIceFrontM2dDakota_25/runme.m	(revision 2110)
@@ -60,5 +60,5 @@
 		%initialize model
 		md=model;
-		md=mesh(md,'DomainOutline.exp',100000);
+		md=mesh(md,'DomainOutline.exp',60000);
 		md=geography(md,'','');
 		md=parameterize(md,'Square.par');
@@ -75,4 +75,6 @@
 		md=tres(md,'dakota');
 		md.results.dakota.importancefactors=importancefactors(md,'drag','max_vel');
+		varargout{1}=md;
+		return;
 
 		%compute fields to be checked