Index: issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp
===================================================================
--- issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp	(revision 3906)
@@ -0,0 +1,79 @@
+/*!\file NodeConnectivityx
+ * \brief: compute node connectivity table, using elements connectivity table.
+ *
+ * For each node, we want to know how many elements are connected to this element, and which they are. 
+ * Given that the 2d meshes we create in ISSM are triangular for now, and they are delaunay conforming, 
+ * each triangle has a minimum angle of 30 degrees, which implies a connectivity <=6. We therefore return 
+ * a nods x 7 connectivity table, with the 7'th column giving us the number of elements connected to each 
+ * row node, and the first 6 columns giving us the elements numbers. 
+ * Amend that: sounds like some triangles get up to 9 connectivity. Take 10 to be on the safe side.
+ * In order to be compatible with matlab output, the connectivity table is given in matlab indexing (starts at 1).
+ */
+
+#include "./NodeConnectivityx.h"
+
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../toolkits/toolkits.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+void	NodeConnectivityx( double** pconnectivity, int* pwidth, double* elements, int nel, int nods){
+
+	int i,j,n;
+	const int maxels=25;
+	const int width=maxels+1;
+
+	/*intermediary: */
+	int     node;
+	int     index;
+	int     num_elements;
+	int     already_plugged=0;
+	double  element;
+
+	/*output: */
+	double* connectivity=NULL;
+
+
+
+	/*Allocate connectivity: */
+	connectivity=(double*)xcalloc(nods*width,sizeof(double));
+
+	/*Go through all elements, and for each elements, plug into the connectivity, all the nodes. 
+	 * If nodes are already plugged into the connectivity, skip them.: */
+	for(n=0;n<nel;n++){
+
+		element=(double)(n+1); //matlab indexing
+
+		for(i=0;i<3;i++){
+		
+			node=(int)*(elements+n*3+i); //already matlab indexed, elements comes directly from the workspace.
+			index=node-1;
+
+			num_elements=(int)*(connectivity+width*index+maxels); //retrieve number of elements already  plugged into the connectivity of this node.
+			
+			already_plugged=0;
+			for(j=0;j<num_elements;j++){
+				if (element==*(connectivity+width*index+j)){
+					already_plugged=1;
+					break;
+				}
+			}
+			if(already_plugged)break;
+
+			/*this elements is not yet plugged  into the connectivity for this node, do it, and increase counter: */
+			*(connectivity+width*index+num_elements)=element;
+			*(connectivity+width*index+maxels)=(double)(num_elements+1);
+			
+		}
+	}
+
+	/*Last check: is the number of elements on last column of the connectivity superior to maxels? If so, then error out and 
+	 * warn the user to increase the connectivity width: */
+	for(i=0;i<nods;i++){
+		if (*(connectivity+width*i+maxels)>maxels)ISSMERROR("%s%g%s"," max connectivity width reached (",*(connectivity+width*i+maxels),")! increase width of connectivity table");
+	}
+
+	/*Assign output pointers: */
+	*pconnectivity=connectivity;
+	*pwidth=width;
+}
Index: issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h
===================================================================
--- issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h	(revision 3906)
+++ issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h	(revision 3906)
@@ -0,0 +1,12 @@
+/*!\file:  NodeConnectivityx.h
+ * \brief header file for node connectivity computation
+ */ 
+
+#ifndef _NODECONNECTIVITYX_H
+#define _NODECONNECTIVITYX_H
+
+/* local prototypes: */
+void	NodeConnectivityx( double** pconnectivity, int* pwidth,double* elements, int nel, int nods);
+
+#endif  /* _NODECONNECTIVITYX_H */
+
Index: issm/trunk/src/c/modules/NormalizeConstraintsx/NormalizeConstraintsx.cpp
===================================================================
--- issm/trunk/src/c/modules/NormalizeConstraintsx/NormalizeConstraintsx.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/NormalizeConstraintsx/NormalizeConstraintsx.cpp	(revision 3906)
@@ -0,0 +1,66 @@
+/*!\file NormalizeConstraintsx
+ * \brief normalize and reduce constraints matrix Rmg to Gmn:
+ */
+
+#include "./NormalizeConstraintsx.h"
+
+#include "../shared/shared.h"
+
+void NormalizeConstraintsx( Mat* pGmn, Mat Rmg, NodeSets* nodesets){
+
+	/*output: */
+	Mat Gmn=NULL;
+
+	/*intermediary: */
+	int i;
+	Mat Rmm=NULL;
+	Mat Rmn=NULL;
+	Mat InvRmm=NULL;
+	int msize,nsize;
+	double* pv_m=NULL;
+	double* row_m=NULL;
+	double* pv_n=NULL;
+	int     null=0;
+
+	if(nodesets){
+		/*Recover data: */
+		msize=nodesets->GetMSize();
+		nsize=nodesets->GetNSize();
+		pv_m=nodesets->GetPV_M();
+		pv_n=nodesets->GetPV_N();
+
+
+		if (msize){
+
+			/*Build row partitioning vector for Rmm and Rmn: */
+			row_m=(double*)xmalloc(msize*sizeof(double));
+			for(i=0;i<msize;i++)row_m[i]=i+1; //matlab indexing for nodesets
+			
+			/*Partition Rmg into Rmm and Rmn: */
+			MatPartition(&Rmm,Rmg,row_m,msize,pv_m,msize); //equivalent of Rmm=Rmg(:,mset);
+			MatPartition(&Rmn,Rmg,row_m,msize,pv_n,nsize);
+			
+			MatView(Rmm,PETSC_VIEWER_STDOUT_WORLD);
+
+			/*Invert Rmm: */
+			MatInvert(&InvRmm,Rmm);
+
+			/*Do Gmn=-(Rmm-1)*Rmn :*/
+			MatScale(InvRmm,-1.0);
+			MatMatMult(InvRmm,Rmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Gmn);
+
+		}
+		else{
+			Gmn=NULL;
+		}
+
+		/*Free ressources:*/
+		xfree((void**)&row_m);
+		MatFree(&Rmm);
+		MatFree(&Rmn);
+		MatFree(&InvRmm);
+	}
+	
+	/*Assign output pointers:*/
+	*pGmn=Gmn;
+}
Index: issm/trunk/src/c/modules/NormalizeConstraintsx/NormalizeConstraintsx.h
===================================================================
--- issm/trunk/src/c/modules/NormalizeConstraintsx/NormalizeConstraintsx.h	(revision 3906)
+++ issm/trunk/src/c/modules/NormalizeConstraintsx/NormalizeConstraintsx.h	(revision 3906)
@@ -0,0 +1,13 @@
+/*!\file:  NormalizeConstraintsx.h
+ * \brief normalize and reduce constraints matrix Rmg to Gmn:
+ */ 
+
+#ifndef _NORMALIZECONSTRAINTSX_H_
+#define _NORMALIZECONSTRAINTSX_H_
+
+#include "../objects/objects.h"
+
+/* local prototypes: */
+void	NormalizeConstraintsx( Mat* pGmn, Mat Rmg, NodeSets* nodesets);
+
+#endif  /* _NORMALIZECONSTRAINTSX_H_ */
Index: issm/trunk/src/c/modules/Orthx/Orthx.cpp
===================================================================
--- issm/trunk/src/c/modules/Orthx/Orthx.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/Orthx/Orthx.cpp	(revision 3906)
@@ -0,0 +1,36 @@
+/*!\file Orthx
+ * \brief orthogonalize searching directions for inverse control methods
+ */
+
+#include "./Orthx.h"
+
+void	Orthx( Vec* pnewgradj, Vec gradj, Vec oldgradj){
+
+	/*output: */
+	Vec newgradj=NULL;
+
+	/*intermediary:*/
+	double norm_new,norm_old,dot_product;;
+
+	/*Initialize output*/
+	VecDuplicate(gradj,&newgradj);
+	VecCopy(gradj,newgradj);
+
+	/*rough orthagonalization
+	gradient=gradient-(gradient'*oldgradient)*oldgradient /norm(oldgradient)^2; */
+	if(oldgradj){
+		VecNorm(oldgradj,NORM_2,&norm_old);
+		VecDot(newgradj,oldgradj,&dot_product);
+		VecAXPY(newgradj, -dot_product/pow(norm_old,2), oldgradj);
+	}
+
+	/*scale to 1: gradient=gradient/max(abs(gradient))*/
+	VecNorm(newgradj,NORM_INFINITY,&norm_new);
+	if (norm_new<=0){
+		ISSMERROR("||∂J/∂α||∞ = 0  gradient is zero");
+	}
+	VecScale(newgradj,1.0/norm_new);
+
+	/*Assign correct pointer*/
+	*pnewgradj=newgradj;
+}
Index: issm/trunk/src/c/modules/Orthx/Orthx.h
===================================================================
--- issm/trunk/src/c/modules/Orthx/Orthx.h	(revision 3906)
+++ issm/trunk/src/c/modules/Orthx/Orthx.h	(revision 3906)
@@ -0,0 +1,16 @@
+/*!\file:  Orthx.h
+ * \brief orthogonalize searching directions for inverse control methods
+ */ 
+
+#ifndef _ORTHX_H
+#define _ORTHX_H
+
+#include "../objects/objects.h"
+#include "../include/include.h"
+#include "../shared/shared.h"
+
+/* local prototypes: */
+void	Orthx( Vec* pnewgradj, Vec gradj, Vec oldgradj);
+
+#endif  /* _ORTHX_H */
+
Index: issm/trunk/src/c/modules/OutputRiftsx/OutputRiftsx.cpp
===================================================================
--- issm/trunk/src/c/modules/OutputRiftsx/OutputRiftsx.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/OutputRiftsx/OutputRiftsx.cpp	(revision 3906)
@@ -0,0 +1,30 @@
+/*!\file OutputRiftsx
+ * \brief: output results from diagnostic solution, for rifts. Notably: fraction of 
+ * melange, and penetration.
+ */
+
+#include "./OutputRiftsx.h"
+
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../toolkits/toolkits.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+void OutputRiftsx( Vec* priftproperties, DataSet* loads, int numrifts){
+
+	/*output: */
+	Vec riftproperties=NULL;
+
+	/*Allocate grad_g: */
+	riftproperties=NewVec(numrifts);
+
+	/*Compute rift properties : */
+	loads->OutputRifts(riftproperties);
+
+	/*Assemble vector: */
+	VecAssemblyBegin(riftproperties);
+	VecAssemblyEnd(riftproperties);
+
+	/*Assign output pointers: */
+	*priftproperties=riftproperties;
+}
Index: issm/trunk/src/c/modules/OutputRiftsx/OutputRiftsx.h
===================================================================
--- issm/trunk/src/c/modules/OutputRiftsx/OutputRiftsx.h	(revision 3906)
+++ issm/trunk/src/c/modules/OutputRiftsx/OutputRiftsx.h	(revision 3906)
@@ -0,0 +1,14 @@
+/*!\file:  OutputRiftsx.h
+ * \brief header file for rift results output.
+ */ 
+
+#ifndef _OUTPUTRIFTSX_H
+#define _OUTPUTRIFTSX_H
+
+#include "../DataSet/DataSet.h"
+
+/* local prototypes: */
+void OutputRiftsx( Vec* priftproperties, DataSet* loads, int numrifts);
+
+#endif  /* _OUTPUTRIFTSX_H */
+
Index: issm/trunk/src/c/modules/PenaltyConstraintsx/PenaltyConstraintsx.cpp
===================================================================
--- issm/trunk/src/c/modules/PenaltyConstraintsx/PenaltyConstraintsx.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/PenaltyConstraintsx/PenaltyConstraintsx.cpp	(revision 3906)
@@ -0,0 +1,50 @@
+/*!\file PenaltyConstraintsx
+ * \brief: set up penalty constraints on loads
+ */
+
+#include "./PenaltyConstraintsx.h"
+#include "./RiftConstraints.h"
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../toolkits/toolkits.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+void PenaltyConstraintsx(int* pconverged, int* pnum_unstable_constraints, DataSet* elements,DataSet* nodes, DataSet* vertices,
+		DataSet* loads,DataSet* materials,  Parameters* parameters,int analysis_type,int sub_analysis_type){
+
+	int i;
+
+	extern int num_procs;
+	extern int my_rank;
+	
+	/*output: */
+	int converged=0;
+	int num_unstable_constraints=0;
+	int min_mechanical_constraints=0;
+
+	/*recover parameters: */
+	parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum);
+
+	/*First, get nodes and loads configured: */
+	elements->Configure(elements, loads, nodes,vertices, materials,parameters);
+	nodes->Configure(elements, loads, nodes,vertices, materials,parameters);
+	loads->Configure(elements, loads, nodes,vertices, materials,parameters);
+
+	/*Do we have penalties linked to rifts? In this case, run our special rifts penalty 
+	 * management routine, otherwise, skip : */
+	if (RiftIsPresent(loads)){
+		RiftConstraints(&converged,&num_unstable_constraints,loads,min_mechanical_constraints,analysis_type,sub_analysis_type);
+	}
+	else if(loads->MeltingIsPresent()){
+		loads->MeltingConstraints(&converged,&num_unstable_constraints,analysis_type,sub_analysis_type);
+	}
+	else{
+		/*Do nothing, no constraints management!:*/
+		num_unstable_constraints=0;
+		converged=1;
+	}
+		
+	/*Assign output pointers: */
+	*pconverged=converged;
+	*pnum_unstable_constraints=num_unstable_constraints;
+}
Index: issm/trunk/src/c/modules/PenaltyConstraintsx/PenaltyConstraintsx.h
===================================================================
--- issm/trunk/src/c/modules/PenaltyConstraintsx/PenaltyConstraintsx.h	(revision 3906)
+++ issm/trunk/src/c/modules/PenaltyConstraintsx/PenaltyConstraintsx.h	(revision 3906)
@@ -0,0 +1,16 @@
+/*!\file:  PenaltyConstraintsx.h
+ * \brief header file for penalty constraints module
+ */ 
+
+#ifndef _PENALTYCONSTRAINTSX_H
+#define _PENALTYCONSTRAINTSX_H
+
+#include "../DataSet/DataSet.h"
+#include "../objects/objects.h"
+
+/* local prototypes: */
+void PenaltyConstraintsx(int* pconverged, int* pnum_unstable_constraints, DataSet* elements,DataSet* nodes, DataSet* vertices,
+		DataSet* loads,DataSet* materials,  Parameters* parameters,int analysis_type,int sub_analysis_type); 
+
+#endif  /* _PENALTYCONSTRAINTSX_H */
+
Index: issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.cpp
===================================================================
--- issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.cpp	(revision 3906)
@@ -0,0 +1,345 @@
+/*!\file RiftConstraints.cpp
+ * \brief: manage penalties for rifts 
+ */
+
+#include "./RiftConstraints.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../include/include.h"
+#include "../shared/shared.h"
+
+#define _ZIGZAGCOUNTER_
+
+int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,int min_mechanical_constraints,int analysis_type,int sub_analysis_type){
+
+	int num_unstable_constraints=0;
+	int converged=0;
+	int potential;
+	extern int my_rank;
+
+	Constrain(&num_unstable_constraints,loads,analysis_type);
+	if(num_unstable_constraints==0)converged=1;
+	
+	
+	if(IsFrozen(loads)){
+		converged=1;
+		num_unstable_constraints=0;
+	}
+	else if(num_unstable_constraints<=min_mechanical_constraints){
+		_printf_("   freezing constraints\n");
+		FreezeConstraints(loads,analysis_type);
+	}
+
+	/*Assign output pointers: */
+	*pconverged=converged;
+	*pnum_unstable_constraints=num_unstable_constraints;
+}
+
+int IsMaterialStable(DataSet* loads,int analysis_type){
+
+	int i;
+	
+	Riftfront* riftfront=NULL;
+	int found=0;
+	int mpi_found=0;
+
+	/*go though loads, and if non-linearity of the material has converged, let all penalties know: */
+	for (i=0;i<loads->Size();i++){
+
+		if(RiftfrontEnum==loads->GetEnum(i)){
+
+			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
+
+			if (riftfront->IsMaterialStable(analysis_type)){
+				found=1;
+				/*do not break! all penalties should get informed the non-linearity converged!*/
+			}
+		}
+	}
+
+	#ifdef _PARALLEL_
+	MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
+	MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);                
+	found=mpi_found;
+	#endif
+
+	return found;
+}
+
+int RiftIsPresent(DataSet* loads){
+
+
+	int i;
+	
+	int found=0;
+	int mpi_found=0;
+
+	/*go though loads, and figure out if one of the loads is a Riftfront: */
+	for (i=0;i<loads->Size();i++){
+
+		if(RiftfrontEnum==loads->GetEnum(i)){
+			found=1;
+			break;
+		}
+	}
+
+	#ifdef _PARALLEL_
+	MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
+	MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);                
+	found=mpi_found;
+	#endif
+
+	return found;
+}
+
+int IsPreStable(DataSet* loads){
+
+
+	int i;
+	
+	Riftfront* riftfront=NULL;
+	int found=0;
+	int mpi_found=0;
+
+	/*go though loads, and figure out if one of the penpair loads is still not stable: */
+	for (i=0;i<loads->Size();i++){
+
+		if(RiftfrontEnum==loads->GetEnum(i)){
+
+			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
+
+			if (riftfront->PreStable()==0){
+				found=1;
+				break;
+			}
+		}
+	}
+
+	#ifdef _PARALLEL_
+	MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
+	MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD);                
+	found=mpi_found;
+	#endif
+
+	if (found){
+		/*We found an unstable constraint. : */
+		return 0;
+	}
+	else{
+		return 1;
+	}
+}
+
+int SetPreStable(DataSet* loads){
+
+
+	int i;
+	
+	Riftfront* riftfront=NULL;
+	int found=0;
+	int mpi_found=0;
+
+	/*go though loads, and set loads to pre stable.:*/
+	for (i=0;i<loads->Size();i++){
+
+		if(RiftfrontEnum==loads->GetEnum(i)){
+
+			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
+			riftfront->SetPreStable();
+		}
+	}
+}
+
+int PreConstrain(int* pnum_unstable_constraints,DataSet* loads,int analysis_type){
+
+	int			i;
+	
+	/* generic object pointer: */
+	Riftfront* riftfront=NULL;
+
+	int unstable;
+	int sum_num_unstable_constraints;
+	int num_unstable_constraints=0;	
+		
+	/*Enforce constraints: */
+	for (i=0;i<loads->Size();i++){
+
+		if (RiftfrontEnum==loads->GetEnum(i)){
+
+			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
+
+			riftfront->PreConstrain(&unstable,analysis_type);
+
+			num_unstable_constraints+=unstable;
+		}
+	}
+
+	#ifdef _PARALLEL_
+	MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
+	MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD);                
+	num_unstable_constraints=sum_num_unstable_constraints;
+	#endif
+	
+	/*Assign output pointers: */
+	*pnum_unstable_constraints=num_unstable_constraints;
+
+}
+
+int Constrain(int* pnum_unstable_constraints,DataSet* loads,int analysis_type){
+
+	int			i;
+	
+	/* generic object pointer: */
+	Riftfront* riftfront=NULL;
+
+	int unstable;
+	int sum_num_unstable_constraints;
+	int num_unstable_constraints=0;	
+		
+	/*Enforce constraints: */
+	for (i=0;i<loads->Size();i++){
+
+		if (RiftfrontEnum==loads->GetEnum(i)){
+
+			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
+
+			riftfront->Constrain(&unstable,analysis_type);
+
+			num_unstable_constraints+=unstable;
+		}
+	}
+
+	#ifdef _PARALLEL_
+	MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
+	MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD);                
+	num_unstable_constraints=sum_num_unstable_constraints;
+	#endif
+	
+	/*Assign output pointers: */
+	*pnum_unstable_constraints=num_unstable_constraints;
+
+}
+
+void FreezeConstraints(DataSet* loads,int analysis_type){
+
+	int			i;
+	
+	/* generic object pointer: */
+	Riftfront* riftfront=NULL;
+
+	/*Enforce constraints: */
+	for (i=0;i<loads->Size();i++){
+
+		if (RiftfrontEnum==loads->GetEnum(i)){
+
+			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
+
+			riftfront->FreezeConstraints(analysis_type);
+
+		}
+	}
+
+}
+
+int IsFrozen(DataSet* loads){
+
+	int			i;
+	
+	/* generic object pointer: */
+	Riftfront* riftfront=NULL;
+	int found=0;
+	int mpi_found=0;
+
+	/*Enforce constraints: */
+	for (i=0;i<loads->Size();i++){
+
+		if (RiftfrontEnum==loads->GetEnum(i)){
+
+			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
+			if (riftfront->IsFrozen()){
+				found=1;
+				break;
+			}
+		}
+	}
+	
+	/*Is there just one found? that would mean we have frozen! : */
+	#ifdef _PARALLEL_
+	MPI_Reduce (&found,&mpi_found,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD );
+	MPI_Bcast(&mpi_found,1,MPI_DOUBLE,0,MPI_COMM_WORLD);                
+	found=mpi_found;
+	#endif
+
+	return found;
+}
+
+int MaxPenetrationInInputs(DataSet* loads,int analysis_type){
+
+	int			i;
+	
+	/* generic object pointer: */
+	Riftfront* riftfront=NULL;
+
+	/*rift penetration: */
+	double max_penetration=0;
+	double mpi_max_penetration;
+	double penetration;
+
+	/*Ok, we are going to find the grid pairs which are not penetrating, even though they 
+	 * are penalised. We will release only the one with has least <0 penetration. : */
+
+	max_penetration=0;
+	for (i=0;i<loads->Size();i++){
+
+		if (RiftfrontEnum==loads->GetEnum(i)){
+
+			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
+
+			riftfront->MaxPenetration(&penetration,analysis_type);
+
+			if (penetration>max_penetration)max_penetration=penetration;
+		}
+	}
+
+	#ifdef _PARALLEL_
+	MPI_Reduce (&max_penetration,&mpi_max_penetration,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD );
+	MPI_Bcast(&mpi_max_penetration,1,MPI_DOUBLE,0,MPI_COMM_WORLD);                
+	max_penetration=mpi_max_penetration;
+	#endif
+
+	/*feed max_penetration to inputs: */
+	loads->UpdateInputsFromVector(&max_penetration,MaxPenetrationEnum,ConstantEnum);
+}
+
+int PotentialUnstableConstraints(DataSet* loads,int analysis_type){
+
+	int			i;
+	
+	/* generic object pointer: */
+	Riftfront* riftfront=NULL;
+
+	/*Ok, we are going to find the grid pairs which are not penetrating, even though they 
+	 * are penalised. We will release only the one with has least <0 penetration. : */
+	int unstable=0;
+	int sum_num_unstable_constraints=0;
+	int num_unstable_constraints=0;
+
+	for (i=0;i<loads->Size();i++){
+
+		if (RiftfrontEnum==loads->GetEnum(i)){
+
+			riftfront=(Riftfront*)loads->GetObjectByOffset(i);
+
+			riftfront->PotentialUnstableConstraint(&unstable,analysis_type);
+
+			num_unstable_constraints+=unstable;
+		}
+	}
+
+	#ifdef _PARALLEL_
+	MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD );
+	MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD);                
+	num_unstable_constraints=sum_num_unstable_constraints;
+	#endif
+
+	return num_unstable_constraints;
+}
Index: issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.h
===================================================================
--- issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.h	(revision 3906)
+++ issm/trunk/src/c/modules/PenaltyConstraintsx/RiftConstraints.h	(revision 3906)
@@ -0,0 +1,33 @@
+/*!\file RiftConstraints.h
+ * \brief: manage penalties for rifts 
+ */
+
+
+#ifndef _RIFTCONSTRAINTS_H_
+#define _RIFTCONSTRAINTS_H_
+
+#include "../objects/objects.h"
+#include "../DataSet/DataSet.h"
+
+int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,int min_mechanical_constraints,int analysis_type,int sub_analysis_type);
+
+int RiftIsPresent(DataSet* loads);
+
+int IsPreStable(DataSet* loads);
+
+int SetPreStable(DataSet* loads);
+
+int PreConstrain(int* pnum_unstable_constraints,DataSet* loads,int analysis_type_enum);
+
+int Constrain(int* pnum_unstable_constraints,DataSet* loads,int analysis_type_enum);
+
+void FreezeConstraints(DataSet* loads,int analysis_type);
+
+int MaxPenetrationInInputs(DataSet* loads,int analysis_type_enum);
+
+int PotentialUnstableConstraints(DataSet* loads,int analysis_type_enum);
+
+int IsMaterialStable(DataSet* loads,int analysis_type_enum);
+
+int IsFrozen(DataSet* loads);
+#endif
Index: issm/trunk/src/c/modules/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp
===================================================================
--- issm/trunk/src/c/modules/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp	(revision 3906)
@@ -0,0 +1,48 @@
+/*!\file PenaltySystemMatricesx
+ * \brief: add penalties to system matrices
+ */
+
+#include "./PenaltySystemMatricesx.h"
+
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../toolkits/toolkits.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+void PenaltySystemMatricesx(Mat Kgg, Vec pg,double* pkmax,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads,DataSet* materials, Parameters* parameters,
+		int kflag,int pflag,int analysis_type,int sub_analysis_type){
+	
+	int i;
+
+	extern int num_procs;
+	extern int my_rank;
+	double kmax;
+	
+	/*First, get elements and loads configured: */
+	elements->Configure(elements,loads, nodes,vertices, materials,parameters);
+	nodes->Configure(elements,loads, nodes,vertices, materials,parameters);
+	loads->Configure(elements, loads, nodes,vertices, materials,parameters);
+	parameters->Configure(elements,loads, nodes,vertices, materials,parameters);
+
+	/*Now, figure out maximum value of K_gg, so that we can penalize it correctly: */
+	MatNorm(Kgg,NORM_INFINITY,&kmax);
+
+	/*Add penalties to stiffnesses, from loads: */
+	if(kflag)loads->PenaltyCreateKMatrix(Kgg,kmax,analysis_type,sub_analysis_type);
+	if(pflag)loads->PenaltyCreatePVector(pg,kmax,analysis_type,sub_analysis_type);
+	
+	/*Assemble matrices: */
+	if(kflag){
+		MatAssemblyBegin(Kgg,MAT_FINAL_ASSEMBLY);
+		MatAssemblyEnd(Kgg,MAT_FINAL_ASSEMBLY);
+		MatCompress(Kgg);
+	}
+	if(pflag){
+		VecAssemblyBegin(pg);
+		VecAssemblyEnd(pg);
+	}
+
+	/*Assign output pointers:*/
+	if(pkmax)*pkmax=kmax;
+
+}
Index: issm/trunk/src/c/modules/PenaltySystemMatricesx/PenaltySystemMatricesx.h
===================================================================
--- issm/trunk/src/c/modules/PenaltySystemMatricesx/PenaltySystemMatricesx.h	(revision 3906)
+++ issm/trunk/src/c/modules/PenaltySystemMatricesx/PenaltySystemMatricesx.h	(revision 3906)
@@ -0,0 +1,16 @@
+/*!\file:  PenaltySystemMatricesx.h
+ * \brief 
+ */ 
+
+#ifndef _PENALTYSYSTEMMATRICESX_H
+#define _PENALTYSYSTEMMATRICESX_H
+
+#include "../DataSet/DataSet.h"
+#include "../objects/objects.h"
+
+/* local prototypes: */
+void PenaltySystemMatricesx(Mat Kgg, Vec pg,double* pkmax, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads,DataSet* materials, Parameters* parameters,
+		int kflag,int pflag,int analysis_type,int sub_analysis_type); 
+
+#endif  /* _PENALTYSYSTEMMATRICESX_H */
+
Index: issm/trunk/src/c/modules/ProcessParamsx/ProcessParamsx.cpp
===================================================================
--- issm/trunk/src/c/modules/ProcessParamsx/ProcessParamsx.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/ProcessParamsx/ProcessParamsx.cpp	(revision 3906)
@@ -0,0 +1,41 @@
+/*!\file ProcessParamsx
+ * \brief: process parameters using partitioning vector. 
+ * Go through all parameters in the 'parameters' dataset. For each parameter that holds a doublevec (ie, a double* vector synchronized across 
+ * the MPI ring of a cluster), partition the vector so that the new node partitioning decided in ModelProcessor is applied. Otherwise, 
+ * parameters coming directly from Matlab would be serially partitioned, which means could not be recognized by each individual node or element. 
+ * The partition needs to be the parallel partitionting.
+ */
+
+#include "./ProcessParamsx.h"
+
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../toolkits/toolkits.h"
+#include "../DataSet/DataSet.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+void ProcessParamsx( Parameters* parameters, Vec  part){
+
+	
+	int     i;
+	double *partition        = NULL;
+	int     numberofvertices;
+	Param  *param            = NULL;
+
+	/*Need number of vertices to repartition DoubleVecParam objects: */
+	parameters->FindParam(&numberofvertices,NumberOfVerticesEnum);
+
+	/*serialize partition vector: */
+	if(part)VecToMPISerial(&partition,part);
+
+	for(i=0;i<parameters->Size();i++){
+
+		param=(Param*)parameters->GetObjectByOffset(i);
+		param->Process(partition,numberofvertices);
+
+	}
+
+	/*Free ressources:*/
+	xfree((void**)&partition);
+
+}
Index: issm/trunk/src/c/modules/ProcessParamsx/ProcessParamsx.h
===================================================================
--- issm/trunk/src/c/modules/ProcessParamsx/ProcessParamsx.h	(revision 3906)
+++ issm/trunk/src/c/modules/ProcessParamsx/ProcessParamsx.h	(revision 3906)
@@ -0,0 +1,14 @@
+/*!\file:  ProcessParamsx.h
+ * \brief header file for processing parameters
+ */ 
+
+#ifndef _PROCESSPARAMSX_H
+#define _PROCESSPARAMSX_H
+
+#include "../DataSet/DataSet.h"
+
+/* local prototypes: */
+void		ProcessParamsx( Parameters* parameters, Vec partition);
+
+#endif  /* _PROCESSPARAMSX_H */
+
Index: issm/trunk/src/c/modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.cpp
===================================================================
--- issm/trunk/src/c/modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.cpp	(revision 3906)
@@ -0,0 +1,38 @@
+/*!\file PropagateFlagsFromConnectivityx
+ */
+
+#include "./PropagateFlagsFromConnectivityx.h"
+
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../toolkits/toolkits.h"
+
+void RecursivePropagation(double* pool, double* connectivity,int index, double* flags);
+
+void PropagateFlagsFromConnectivityx( double* pool, double* connectivity,int index, double* flags){
+
+	/*Call recursive propagation routine: */
+	RecursivePropagation(pool, connectivity,index, flags);
+}
+
+
+void RecursivePropagation(double* pool, double* connectivity, int index, double* flags){
+
+	int i;
+	int newel;
+
+	/*if this element (index) belongs to the pool already, skip: */
+	if(pool[index-1])return;
+
+	/*if this element does not belong to the flags set, skip: */
+	if(flags[index-1]==0)return;
+
+	/*put this element (index), which belongs to the flags, into the pool: */
+	pool[index-1]=1;
+
+	/*now, propagate recursively using connectivity of this element: */
+	for(i=0;i<3;i++){
+		newel=(int)*(connectivity+(index-1)*3+i);
+		RecursivePropagation(pool, connectivity, newel, flags);
+	}
+}
Index: issm/trunk/src/c/modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h
===================================================================
--- issm/trunk/src/c/modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h	(revision 3906)
+++ issm/trunk/src/c/modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h	(revision 3906)
@@ -0,0 +1,13 @@
+/*!\file:  PropagateFlagsFromConnectivityx.h
+ */ 
+
+#ifndef _PROPAGATEFLAGSFROMCONNECTIVITYX_H
+#define _PROPAGATEFLAGSFROMCONNECTIVITYX_H
+
+#include "../DataSet/DataSet.h"
+
+/* local prototypes: */
+void PropagateFlagsFromConnectivityx( double* pool, double* connectivity,int index, double* flags);
+
+#endif  /* _PROPAGATEFLAGSFROMCONNECTIVITYX_H */
+
Index: issm/trunk/src/c/modules/Qmux/DakotaResponses.cpp
===================================================================
--- issm/trunk/src/c/modules/Qmux/DakotaResponses.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/Qmux/DakotaResponses.cpp	(revision 3906)
@@ -0,0 +1,338 @@
+/*!\file:  DakotaResponses.cpp
+ * \brief  compute dakota responses, using a list of response descriptors.
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "../DataSet/DataSet.h"    
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../MassFluxx/MassFluxx.h"
+#include "../Misfitx/Misfitx.h"
+#include "../modules.h"
+
+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;
+	int found=0;
+	char* response_descriptor=NULL;
+	int numberofnodes;
+	extern int my_rank;
+
+	/*some data needed across the responses: */
+	model->FindParam(&numberofnodes,NumberOfNodesEnum);
+
+
+	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, mass_flux'
+
+		if(strcmp(response_descriptor,"min_vel")==0){
+			double* vel=NULL;
+			double min_vel=0;
+		
+			found=processed_results->FindResult((void*)&vel,"vel");
+			if(!found)ISSMERROR(" could not find vel to compute min_vel");
+
+			min_vel=vel[0];
+			for(j=1;j<numberofnodes;j++){
+				if (vel[j]<min_vel)min_vel=vel[j];
+			}
+
+			if(my_rank==0)responses[i]=min_vel;
+			
+			/*Free ressources:*/
+			xfree((void**)&vel);
+
+			
+		}
+		else if(strcmp(response_descriptor,"max_vel")==0){
+			double* vel=NULL;
+			double max_vel=0;
+
+			found=processed_results->FindResult((void*)&vel,"vel");
+			if(!found)ISSMERROR(" could not find vel to compute max_vel");
+
+			max_vel=vel[0];
+			for(j=1;j<numberofnodes;j++){
+				if (vel[j]>max_vel)max_vel=vel[j];
+			}
+			if(my_rank==0)responses[i]=max_vel;
+			
+			/*Free ressources:*/
+			xfree((void**)&vel);
+
+		}
+		else if(strcmp(response_descriptor,"min_vx")==0){
+			double* vx=NULL;
+			double min_vx=0;
+			
+			found=processed_results->FindResult((void*)&vx,"vx");
+			if(!found)ISSMERROR(" could not find vx to compute min_vx");
+
+			min_vx=vx[0];
+			for(j=1;j<numberofnodes;j++){
+				if (vx[j]<min_vx)min_vx=vx[j];
+			}
+			if(my_rank==0)responses[i]=min_vx;
+			
+			/*Free ressources:*/
+			xfree((void**)&vx);
+
+		}
+		else if(strcmp(response_descriptor,"max_vx")==0){
+			double* vx=NULL;
+			double max_vx=0;
+			
+			found=processed_results->FindResult((void*)&vx,"vx");
+			if(!found)ISSMERROR(" could not find vx to compute max_vx");
+
+			max_vx=vx[0];
+			for(j=1;j<numberofnodes;j++){
+				if (vx[j]>max_vx)max_vx=vx[j];
+			}
+			if(my_rank==0)responses[i]=max_vx;
+			
+			/*Free ressources:*/
+			xfree((void**)&vx);
+
+		}
+		else if(strcmp(response_descriptor,"max_abs_vx")==0){
+			double* vx=NULL;
+			double max_abs_vx=0;
+			
+			found=processed_results->FindResult((void*)&vx,"vx");
+			if(!found)ISSMERROR(" could not find vx to compute max_abs_vx");
+
+			max_abs_vx=fabs(vx[0]);
+			for(j=1;j<numberofnodes;j++){
+				if (fabs(vx[j])>max_abs_vx)max_abs_vx=fabs(vx[j]);
+			}
+			if(my_rank==0)responses[i]=max_abs_vx;
+			
+			/*Free ressources:*/
+			xfree((void**)&vx);
+
+		}
+		else if(strcmp(response_descriptor,"min_vy")==0){
+			double* vy=NULL;
+			double min_vy=0;
+			
+			found=processed_results->FindResult((void*)&vy,"vy");
+			if(!found)ISSMERROR(" could not find vy to compute min_vy");
+
+			min_vy=vy[0];
+			for(j=1;j<numberofnodes;j++){
+				if (vy[j]<min_vy)min_vy=vy[j];
+			}
+			if(my_rank==0)responses[i]=min_vy;
+			
+			/*Free ressources:*/
+			xfree((void**)&vy);
+
+		}
+		else if(strcmp(response_descriptor,"max_vy")==0){
+			double* vy=NULL;
+			double max_vy=0;
+			
+			found=processed_results->FindResult((void*)&vy,"vy");
+			if(!found)ISSMERROR(" could not find vy to compute max_vy");
+
+			max_vy=vy[0];
+			for(j=1;j<numberofnodes;j++){
+				if (vy[j]>max_vy)max_vy=vy[j];
+			}
+			if(my_rank==0)responses[i]=max_vy;
+			
+			/*Free ressources:*/
+			xfree((void**)&vy);
+
+		}
+		else if(strcmp(response_descriptor,"max_abs_vy")==0){
+			double* vy=NULL;
+			double max_abs_vy=0;
+			
+			found=processed_results->FindResult((void*)&vy,"vy");
+			if(!found)ISSMERROR(" could not find vy to compute max_abs_vy");
+
+			max_abs_vy=fabs(vy[0]);
+			for(j=1;j<numberofnodes;j++){
+				if (fabs(vy[j])>max_abs_vy)max_abs_vy=fabs(vy[j]);
+			}
+			if(my_rank==0)responses[i]=max_abs_vy;
+			
+			/*Free ressources:*/
+			xfree((void**)&vy);
+
+		}
+		else if(strcmp(response_descriptor,"min_vz")==0){
+			double* vz=NULL;
+			double min_vz=0;
+			
+			found=processed_results->FindResult((void*)&vz,"vz");
+			if(!found)ISSMERROR(" could not find vz to compute min_vz");
+
+			min_vz=vz[0];
+			for(j=1;j<numberofnodes;j++){
+				if (vz[j]<min_vz)min_vz=vz[j];
+			}
+			if(my_rank==0)responses[i]=min_vz;
+			
+			/*Free ressources:*/
+			xfree((void**)&vz);
+
+		}
+		else if(strcmp(response_descriptor,"max_vz")==0){
+			double* vz=NULL;
+			double max_vz=0;
+			
+			found=processed_results->FindResult((void*)&vz,"vz");
+			if(!found)ISSMERROR(" could not find vz to compute max_vz");
+
+			max_vz=vz[0];
+			for(j=1;j<numberofnodes;j++){
+				if (vz[j]>max_vz)max_vz=vz[j];
+			}
+			if(my_rank==0)responses[i]=max_vz;
+			
+			/*Free ressources:*/
+			xfree((void**)&vz);
+
+		}
+		else if(strcmp(response_descriptor,"max_abs_vz")==0){
+			double* vz=NULL;
+			double max_abs_vz=0;
+			
+			found=processed_results->FindResult((void*)&vz,"vz");
+			if(!found)ISSMERROR(" could not find vz to compute max_abs_vz");
+
+			max_abs_vz=fabs(vz[0]);
+			for(j=1;j<numberofnodes;j++){
+				if (fabs(vz[j])>max_abs_vz)max_abs_vz=fabs(vz[j]);
+			}
+			if(my_rank==0)responses[i]=max_abs_vz;
+			
+			/*Free ressources:*/
+			xfree((void**)&vz);
+		}
+		else if(strcmp(response_descriptor,"misfit")==0){
+			
+			int isstokes,ismacayealpattyn,ishutter;
+			FemModel* femmodel=NULL;
+			double J=0;
+			int numberofdofspernode,numberofnodes;
+			Vec u_g=NULL;
+			double* u_g_double=NULL;
+			double* vx=NULL;
+			double* vy=NULL;
+			double* vz=NULL;
+			double* fit=NULL;
+
+			/*retrieve active fem model: */
+			model->FindParam(&isstokes,IsStokesEnum);
+			model->FindParam(&ismacayealpattyn,IsMacAyealPattynEnum);
+			model->FindParam(&ishutter,IsHutterEnum);
+
+			if(isstokes){
+				femmodel=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum);
+			}
+			if(ismacayealpattyn){
+				femmodel=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum);
+			}
+			if(ishutter){
+				femmodel=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum);
+			}	
+
+
+			/*Recover some parameters: */
+			femmodel->parameters->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum);
+			femmodel->parameters->FindParam(&numberofnodes,NumberOfNodesEnum);
+			femmodel->parameters->FindParam(&fit,NULL,NULL,FitEnum);
+
+			/*Recover velocity: */
+			found=results->FindResult(&u_g,"u_g");
+			VecToMPISerial(&u_g_double,u_g);
+			if(!found)ISSMERROR(" could not find velocity to compute misfit");
+
+			SplitSolutionVectorx(u_g,numberofnodes,numberofdofspernode,&vx,&vy,&vz);
+
+			/*Add to inputs: */
+			femmodel->elements->UpdateInputsFromVector(vx,VxEnum,VertexEnum);
+			femmodel->elements->UpdateInputsFromVector(vy,VyEnum,VertexEnum);
+			femmodel->elements->UpdateInputsFromVector(vz,VzEnum,VertexEnum);
+			femmodel->elements->UpdateInputsFromVector(&fit[0],FitEnum,ConstantEnum);
+
+			/*Compute misfit: */
+			Misfitx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,analysis_type,sub_analysis_type);
+			
+
+			if(my_rank==0)responses[i]=J;
+
+			/*Some cleanup: */
+			VecFree(&u_g);
+			xfree((void**)&u_g_double);
+			xfree((void**)&fit);
+
+		}
+		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)ISSMERROR(" could not find velocity to compute mass_flux");
+			VecToMPISerial(&ug_serial,ug);
+		
+			/*retrieve active fem model: */
+			model->FindParam(&isstokes,IsStokesEnum);
+			model->FindParam(&ismacayealpattyn,IsMacAyealPattynEnum);
+			model->FindParam(&ishutter,IsHutterEnum);
+
+			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: */
+			femmodel->parameters->FindParam(&segments,&num_segments,QmuMassFluxSegmentsEnum);
+
+			/*call mass flux module: */
+			MassFluxx(&mass_flux,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,segments,num_segments,ug_serial);
+			
+			if(my_rank==0)responses[i]=mass_flux;
+			
+			/*Free ressources:*/
+			VecFree(&ug);
+			xfree((void**)&ug_serial);
+			xfree((void**)&segments);
+		}
+		else{
+			if(my_rank==0)printf("%s%s%s"," response descriptor : ",response_descriptor," not supported yet!");
+			ISSMERROR("%s%s%s"," response descriptor : ",response_descriptor," not supported yet!");
+		}
+	}
+
+}
Index: issm/trunk/src/c/modules/Qmux/Qmux.cpp
===================================================================
--- issm/trunk/src/c/modules/Qmux/Qmux.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/Qmux/Qmux.cpp	(revision 3906)
@@ -0,0 +1,145 @@
+/*!\file:  Qmux.cpp
+ * \brief: wrapper to the Dakota capabilities. qmu fires up Dakota, and registers a Dakota Pluggin
+ * which will be in charge of running the solution sequences repeteadly, to garner statistics. 
+ *
+ * This routine deals with running ISSM and Dakota in library mode. In library mode, Dakota does not 
+ * run as an execuatble. Its capabilities are linked into the ISSM software. ISSM calls dakota routines 
+ * directly from the dakota library. qmu.cpp is the code that is in charge of calling those routines. 
+ *
+ * Dakota has its own way of running in parallel (for embarassingly parallel jobs). We do not want that, 
+ * as ISSM knows exactly how to run "really parallel" jobs that use all CPUS. To bypass Dakota's parallelism, 
+ * we overloaded the constructor for the parallel library (see the Dakota patch in the externalpackages/dakota
+ * directory). This overloaded constructor fires up Dakota serially on CPU 0 only! We take care of broadcasting 
+ * to the other CPUS, hence ISSM is running in parallel, and Dakota serially on CPU0. 
+ *
+ * Now, how does CPU 0 drive all other CPUS to carry out sensitivity analysese? By synchronizing its call to 
+ * our ISSM cores (diagnostic_core, thermal_core, transient_core, etc ...) on CPU 0 with all other CPUS. 
+ * This explains the structure of qmu.cpp, where cpu 0 runs Dakota, the Dakota pluggin fires up SpawnCore.cpp, 
+ * while the other CPUS are waiting for a broadcast from CPU0, once they get it, they also fire up 
+ * SpawnCore. In the end, SpawnCore is fired up on all CPUS, with CPU0 having Dakota inputs, that it will 
+ * broacast to other CPUS. 
+ *
+ * Now, how does dakota call the SpawnCore routine? The SpawnCore is embedded into the DakotaPlugin object 
+ * which is derived from the Direct Interface Dakota objct. This is the only way to run Dakota in library 
+ * mode (see their developper guide for more info). Dakota registers the DakotaPlugin object into its own 
+ * database, and calls on the embedded SpawnCore from CPU0. 
+ *
+ */ 
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./Qmux.h"
+
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../toolkits/toolkits.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../objects/DakotaPlugin.h"
+
+#ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in.
+#include "ParallelLibrary.H"
+#include "ProblemDescDB.H"
+#include "DakotaStrategy.H"
+#include "DakotaModel.H"
+#include "DakotaInterface.H"
+
+#endif
+
+#ifdef _SERIAL_
+void Qmux(mxArray* model,int analysis_type,int sub_analysis_type,char* dakota_input_file,char* dakota_output_file,char* dakota_error_file){
+#else
+void Qmux(Model* model,int analysis_type,int sub_analysis_type){
+#endif
+
+
+	#ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in.
+
+	#ifdef _PARALLEL_
+	char* dakota_input_file=NULL;
+	char* dakota_output_file=NULL;
+	char* dakota_error_file=NULL;
+	extern int my_rank;
+	#endif
+	int status=0;
+	Dakota::ModelLIter ml_iter;
+
+	#ifdef _PARALLEL_
+	/*Recover dakota_input_file, dakota_output_file and dakota_error_file, in the parameters dataset in parallel */
+	model->FindParam(&dakota_input_file,QmuInNameEnum);
+	model->FindParam(&dakota_output_file,QmuOutNameEnum);
+	model->FindParam(&dakota_error_file,QmuErrNameEnum);
+	#endif
+
+	#ifdef _PARALLEL_
+	if(my_rank==0){
+	#endif
+	
+		// Instantiate/initialize the parallel library and problem description
+		// database objects.
+		#ifdef _SERIAL_
+			Dakota::ParallelLibrary parallel_lib; //use Dakota's standard library mode constructor
+		#else
+			Dakota::ParallelLibrary parallel_lib("serial"); //use our own ISSM Dakota library mode constructor, which only fires up Dakota on CPU 0. 
+		#endif
+		Dakota::ProblemDescDB problem_db(parallel_lib); 
+
+		// Manage input file parsing, output redirection, and restart processing
+		// without a CommandLineHandler.  This version relies on parsing of an
+		// input file.
+		problem_db.manage_inputs(dakota_input_file);
+		// specify_outputs_restart() is only necessary if specifying non-defaults
+		parallel_lib.specify_outputs_restart(dakota_output_file,dakota_error_file,NULL,NULL);
+
+		// Instantiate the Strategy object (which instantiates all Model and
+		// Iterator objects) using the parsed information in problem_db.
+		Dakota::Strategy selected_strategy(problem_db);
+
+		// convenience function for iterating over models and performing any
+		// interface plug-ins
+		Dakota::ModelList& models = problem_db.model_list();
+
+		for (ml_iter = models.begin(); ml_iter != models.end(); ml_iter++) {
+
+			Dakota::Interface& interface = ml_iter->interface();
+
+			//set DB nodes to the existing Model specification
+			problem_db.set_db_model_nodes(ml_iter->model_id());
+
+			// Serial case: plug in derived Interface object without an analysisComm
+			interface.assign_rep(new SIM::DakotaPlugin(problem_db,(void*)model,analysis_type,sub_analysis_type), false);
+		}
+	
+		// Execute the strategy
+		problem_db.lock(); // prevent run-time DB queries
+		selected_strategy.run_strategy();
+		
+		#ifdef _PARALLEL_
+		//Warn other cpus that we are done running the dakota iterator, by setting the counter to -1:
+		SpawnCore(NULL,0, NULL,NULL,0,model,analysis_type,sub_analysis_type,-1);
+		#endif
+
+	#ifdef _PARALLEL_
+	}
+	else{
+
+		for(;;){
+			if(!SpawnCore(NULL,0, NULL,NULL,0,model,analysis_type,sub_analysis_type,0))break; //counter came in at -1 on cpu0, bail out.
+		}
+	}
+	#endif
+
+	/*Free ressources:*/
+	#ifdef _PARALLEL_
+	xfree((void**)&dakota_input_file);
+	xfree((void**)&dakota_error_file);
+	xfree((void**)&dakota_output_file);
+	#endif
+
+
+	#endif
+
+}
Index: issm/trunk/src/c/modules/Qmux/Qmux.h
===================================================================
--- issm/trunk/src/c/modules/Qmux/Qmux.h	(revision 3906)
+++ issm/trunk/src/c/modules/Qmux/Qmux.h	(revision 3906)
@@ -0,0 +1,24 @@
+/*!\file:  Qmux.h
+ * \brief header file for Qmu engine using Dakota
+ */ 
+
+#ifndef _QMUX_H
+#define _QMUX_H
+
+#include "../DataSet/DataSet.h"
+#include "../objects/objects.h"
+
+/* local prototypes: */
+int SpawnCore(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, void* model,int analysis_type,int sub_analysis_type,int counter);
+#ifdef _SERIAL_
+void Qmux(mxArray* model,int analysis_type,int sub_analysis_type,char* dakota_input_file,char* dakota_output_file,char* dakota_error_file);
+void SpawnCoreSerial(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, mxArray* model,int analysis_type,int sub_analysis_type,int counter);
+#else
+void Qmux(Model* model,int analysis_type,int sub_analysis_type);
+void SpawnCoreParallel(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, Model* model,int analysis_type,int sub_analysis_type,int counter);
+void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Model* model, DataSet* results,DataSet* processed_results,int analysis_type,int sub_analysis_type);
+#endif
+
+
+#endif  /* _QMUX_H */
+
Index: issm/trunk/src/c/modules/Qmux/SpawnCore.cpp
===================================================================
--- issm/trunk/src/c/modules/Qmux/SpawnCore.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/Qmux/SpawnCore.cpp	(revision 3906)
@@ -0,0 +1,32 @@
+/*!\file:  SpawnCore.cpp
+ * \brief: branch into SpawnCoreMatlab and SpawnCoreParallel.
+ */
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+
+#include "../objects/objects.h"
+#include "../io/io.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../shared/shared.h"
+#include "./Qmux.h"
+#include "../include/include.h"
+
+int SpawnCore(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, void* model,int analysis_type,int sub_analysis_type,int counter){
+
+	/*Branch into a serial SpawnCore and a parallel SpawnCore: */
+	#ifdef _SERIAL_
+	SpawnCoreSerial(responses, numresponses, variables, variables_descriptors,numvariables, (mxArray*)model, analysis_type,sub_analysis_type,counter);
+	#else
+	/*Call SpawnCoreParallel unless counter=-1 on cpu0, in which case, bail out and return 0: */
+	MPI_Bcast(&counter,1,MPI_INT,0,MPI_COMM_WORLD); if(counter==-1)return 0;
+	
+	SpawnCoreParallel(responses, numresponses, variables, variables_descriptors,numvariables, (Model*)model, analysis_type,sub_analysis_type,counter);
+	#endif
+
+	return 1;
+}
Index: issm/trunk/src/c/modules/Qmux/SpawnCoreParallel.cpp
===================================================================
--- issm/trunk/src/c/modules/Qmux/SpawnCoreParallel.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/Qmux/SpawnCoreParallel.cpp	(revision 3906)
@@ -0,0 +1,158 @@
+/*!\file:  SpawnCoreParallel.cpp
+ * \brief: run core ISSM solution using Dakota inputs coming from CPU 0.
+ * \sa qmu.cpp DakotaPlugin.cpp
+ *
+ * This routine needs to be understood simultaneously with qmu.cpp and DakotaPlugin. 
+ * SpawnCoreParallel is called by all CPUS, with CPU 0 holding Dakota variable values, along 
+ * with variable descriptors. 
+ *
+ * SpawnCoreParallel takes care of broadcasting the variables and their descriptors across the MPI 
+ * ring. Once this is done, we use the variables to modify the inputs for the solution core. 
+ * For ex, if "rho_ice" is provided, for ex 920, we include "rho_ice" in the inputs, then 
+ * call the core with the modified inputs. This is the way we get Dakota to explore the parameter 
+ * spce of the core. 
+ *
+ * Once the core is called, we process the results of the core, and using the processed results, 
+ * we compute response functions. The responses are computed on all CPUS, but they are targeted 
+ * for CPU 0, which will get these values back to the Dakota engine. 
+ *
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+
+#include "../objects/objects.h"
+#include "../io/io.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../shared/shared.h"
+#include "./Qmux.h"
+#include "../include/include.h"
+#include "../solutions/solutions.h"
+
+void SpawnCoreParallel(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, Model* model,int analysis_type,int sub_analysis_type,int counter){
+
+	int i;
+	
+	/*output from core solutions: */
+	DataSet* results=NULL;
+	DataSet* processed_results=NULL;
+
+	char** responses_descriptors=NULL;
+	int    num_responses_descriptors;
+	char*  string=NULL;
+	int    string_length;
+	double* qmu_part=NULL;
+	int     qmu_npart;
+	int     verbose=0;
+	int     dummy;
+
+	extern int my_rank;
+	
+	
+	/*synchronize all cpus, as CPU 0 is probably late (it is starting the entire dakota strategy!) : */
+	MPI_Barrier(MPI_COMM_WORLD);
+	
+	/*some parameters needed: */
+	model->FindParam(&verbose,VerboseEnum);
+		
+	/*First off, recover the response descriptors for the response functions: */
+	model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum)->parameters->FindParam(&responses_descriptors,&num_responses_descriptors,ResponseDescriptorsEnum);
+
+	/*Recover partitioning for dakota: */
+	model->FindParam(&qmu_npart,QmuNPartEnum);
+	model->FindParam(&qmu_part,&dummy,QmuPartEnum);
+
+	/*broadcast variables: only cpu 0 has correct values*/
+	MPI_Bcast(&numvariables,1,MPI_INT,0,MPI_COMM_WORLD); 
+	if(my_rank!=0)variables=(double*)xmalloc(numvariables*sizeof(double));
+	MPI_Bcast(variables,numvariables,MPI_DOUBLE,0,MPI_COMM_WORLD); 
+
+	/*broadcast variables_descriptors: */
+	if(my_rank!=0){
+		variables_descriptors=(char**)xmalloc(numvariables*sizeof(char*));
+	}
+	for(i=0;i<numvariables;i++){
+		if(my_rank==0){
+			string=variables_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)variables_descriptors[i]=string;
+	}
+
+	/*broadcast numresponses: */
+	MPI_Bcast(&numresponses,1,MPI_INT,0,MPI_COMM_WORLD); 
+
+	_printf_("qmu iteration: %i\n",counter);
+
+	/*Modify core inputs in objects contained in model, to reflect the dakota variables inputs: */
+	model->UpdateFromDakota(variables,variables_descriptors,numvariables,model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum)->parameters,qmu_part,qmu_npart); //diagnostic horiz model is the one holding the parameters for Dakota.
+
+	/*Run the analysis core solution sequence: */
+	if(analysis_type==DiagnosticAnalysisEnum){
+			
+		if(verbose)_printf_("Starting diagnostic core\n");
+
+		results=diagnostic_core(model);
+
+	}
+	else if(analysis_type==ThermalAnalysisEnum){
+		
+		if(verbose)_printf_("Starting thermal core\n");
+		results=thermal_core(model);
+
+	}
+	else if(analysis_type==PrognosticAnalysisEnum){
+
+		if(verbose)_printf_("Starting prognostic core\n");
+		results=prognostic_core(model);
+
+	}
+	else if(analysis_type==TransientAnalysisEnum){
+
+		if(verbose)_printf_("Starting transient core\n");
+		results=transient_core(model);
+
+	}
+	else ISSMERROR("%s%i%s%i%s"," analysis_type ",analysis_type," and sub_analysis_type ",sub_analysis_type," not supported yet!");
+	
+		
+
+	/*Now process the outputs, before computing the dakota responses: */
+	if(verbose)_printf_("process results:\n");
+
+	ProcessResults(&processed_results,results,model,analysis_type); 
+
+	/*compute responses on cpu 0: dummy for now! */
+	if(verbose)_printf_("compute dakota responses:\n");
+	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
+	if(my_rank!=0){
+		xfree((void**)&variables);
+		for(i=0;i<numvariables;i++){
+			string=variables_descriptors[i];
+			xfree((void**)&string);
+		}
+		xfree((void**)&variables_descriptors);
+	}
+	//responses descriptors
+	for(i=0;i<numresponses;i++){
+		string=responses_descriptors[i];
+		xfree((void**)&string);
+	}
+	//rest of dynamic allocations.
+	xfree((void**)&responses_descriptors);
+	xfree((void**)&qmu_part);
+}
+
Index: issm/trunk/src/c/modules/Qmux/SpawnCoreSerial.cpp
===================================================================
--- issm/trunk/src/c/modules/Qmux/SpawnCoreSerial.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/Qmux/SpawnCoreSerial.cpp	(revision 3906)
@@ -0,0 +1,79 @@
+/*!\file:  SpawnCoreSerial.cpp
+ * \brief: run core ISSM solution using Dakota inputs. Call the Serial core solution, using mexCallMATLAB
+ * \sa SpawnCore.cpp SpawnCoreParallel.cpp
+ *
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+
+#include "../objects/objects.h"
+#include "../io/io.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../shared/shared.h"
+#include "./Qmux.h"
+#include "../include/include.h"
+
+void SpawnCoreSerial(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, mxArray* model,int analysis_type,int sub_analysis_type,int counter){
+
+	int i;
+	
+	//inputs to matlab routine
+	mxArray* mxvariables=NULL;
+	double*  variables_copy=NULL;
+	mxArray* mxvariabledescriptors=NULL;
+	mxArray* mxanalysis_type=NULL;
+	mxArray* mxsub_analysis_type=NULL;
+	mxArray* mxcounter=NULL;
+	mwSize   dims[2]={0};
+
+	//mexCallMATLAB arrays
+	mxArray* array[7];
+
+	//output from SpawnCore in matlab routine.
+	mxArray* mxresponses=NULL;
+
+	/*Create variables and variable descriptors mxArrays that we will feed to the core solution for update of the inputs: */
+	mxvariables=mxCreateDoubleMatrix(numvariables,1,mxREAL);
+	variables_copy=(double*)xmalloc(numvariables*sizeof(double));
+	memcpy(variables_copy,variables,numvariables*sizeof(double));
+	mxSetPr(mxvariables,variables_copy);
+
+	dims[0]=numvariables;
+	dims[1]=1;
+	mxvariabledescriptors=mxCreateCellArray(2,dims);
+	for(i=0;i<numvariables;i++){
+		mxSetCell(mxvariabledescriptors,i,mxCreateString(variables_descriptors[i]));
+	}
+
+	mxanalysis_type=mxCreateDoubleScalar((double)analysis_type);
+	mxsub_analysis_type=mxCreateDoubleScalar((double)sub_analysis_type);
+	mxcounter=mxCreateDoubleScalar((double)counter);
+
+	//call SpwanCore matlab routine.
+	array[0]=model;
+	array[2]=mxvariables;
+	array[3]=mxvariabledescriptors;
+	array[4]=mxanalysis_type;
+	array[5]=mxsub_analysis_type;
+	array[6]=mxcounter;
+
+	mexCallMATLAB(1,&mxresponses,7,array,"SpawnCore");
+
+	/*copy responses back to dakota: */
+	memcpy(responses,mxGetPr(mxresponses),numresponses*sizeof(double));
+
+	//destroy constructed arrays: 
+	mxDestroyArray(mxvariables);
+	mxDestroyArray(mxvariabledescriptors);
+	mxDestroyArray(mxresponses);
+	mxDestroyArray(mxanalysis_type);
+	mxDestroyArray(mxsub_analysis_type);
+	mxDestroyArray(mxcounter);
+
+}
+
Index: issm/trunk/src/c/modules/Reduceloadfromgtofx/Reduceloadfromgtofx.cpp
===================================================================
--- issm/trunk/src/c/modules/Reduceloadfromgtofx/Reduceloadfromgtofx.cpp	(revision 3906)
+++ issm/trunk/src/c/modules/Reduceloadfromgtofx/Reduceloadfromgtofx.cpp	(revision 3906)
@@ -0,0 +1,93 @@
+/*!\file Reduceloadfromgtofx
+ * \brief reduce loads from g set to f set 
+ */
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./Reduceloadfromgtofx.h"
+
+void	Reduceloadfromgtofx( Vec* ppf, Vec pg, Mat G_mn, Mat Kfs, Vec y_s, NodeSets* nodesets){
+
+	/*output: */
+	Vec pf=NULL;
+
+	/*intermediary*/
+	Vec pn=NULL;
+	Vec pm=NULL;
+	Mat tG_mn=NULL;
+	int tG_mnm,tG_mnn;
+	Vec tG_mnpm=NULL;
+	Vec Kfsy_s=NULL;
+	int Kfsm,Kfsn;
+	PetscScalar a;
+	
+
+	if(!pg){
+		pf=NULL;
+	}
+	else{
+
+		/* Reduce pg to pn:*/
+		if (nodesets->GetMSize()){
+
+			VecPartition(&pn, pg, nodesets->GetPV_N(), nodesets->GetNSize());
+			VecPartition(&pm, pg, nodesets->GetPV_M(), nodesets->GetMSize());
+		
+			/*Create G_mn'*/
+			#if _PETSC_VERSION_ == 2
+			MatTranspose(G_mn,&tG_mn);  
+			#else
+			MatTranspose(G_mn,MAT_INITIAL_MATRIX,&tG_mn);  
+			#endif
+
+			/*Create G_mn' * pm*/
+			MatGetLocalSize(tG_mn,&tG_mnm,&tG_mnn);
+			tG_mnpm=NewVecFromLocalSize(tG_mnm);
+			MatMultPatch(tG_mn,pm,tG_mnpm);
+
+			/*add pn to G_mn' * pm*/ 
+			a=1;
+			VecAXPY(pn,a,tG_mnpm);  
+		}
+		else{
+			VecDuplicate(pg,&pn);  
+			VecCopy(pg,pn);  
+		}
+		
+		/*% Reduce pn to pf:*/
+
+		if (nodesets->GetSSize()){
+
+			VecPartition(&pf, pn, nodesets->GetPV_F(),nodesets->GetFSize());
+
+			/*pf = pf - Kfs * y_s;*/
+			MatGetLocalSize(Kfs,&Kfsm,&Kfsn);
+			Kfsy_s=NewVecFromLocalSize(Kfsm);
+			MatMultPatch(Kfs,y_s,Kfsy_s);
+
+			a=-1;
+			VecAXPY(pf,a,Kfsy_s);  
+		}
+		else{
+			/*pf=pn*/
+			VecDuplicate(pn,&pf);  
+			VecCopy(pn,pf);  
+		}
+	}
+	
+	
+	/*Assign correct pointer*/
+	*ppf=pf;
+
+	/*Free ressources and return*/
+	VecFree(&pn);
+	VecFree(&pm);
+	MatFree(&tG_mn);
+	VecFree(&tG_mnpm);
+	VecFree(&Kfsy_s);
+
+}
Index: issm/trunk/src/c/modules/Reduceloadfromgtofx/Reduceloadfromgtofx.h
===================================================================
--- issm/trunk/src/c/modules/Reduceloadfromgtofx/Reduceloadfromgtofx.h	(revision 3906)
+++ issm/trunk/src/c/modules/Reduceloadfromgtofx/Reduceloadfromgtofx.h	(revision 3906)
@@ -0,0 +1,14 @@
+/*!\file:  Reduceloadfromgtofx.h
+ * \brief reduce load from g set to f set
+ */ 
+
+#ifndef _REDUCELOADFROMGTOFX_H
+#define _REDUCELOADFROMGTOFX_H
+
+#include "../objects/objects.h"
+
+/* local prototypes: */
+void	Reduceloadfromgtofx( Vec* ppf, Vec pg, Mat Gmn, Mat Kfs, Vec ys, NodeSets* nodesets);
+
+#endif  /* _REDUCELOADFROMGTOFX_H */
+
