Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 25999)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 26000)
@@ -641,5 +641,5 @@
 	./modules/NodeConnectivityx/NodeConnectivityx.cpp \
 	./modules/ElementConnectivityx/ElementConnectivityx.cpp \
-	./modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.cpp 
+	./modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.cpp
 
 if CHACO
@@ -783,5 +783,5 @@
 endif
 
-issm_SOURCES = main/issm.cpp 
+issm_SOURCES = main/issm.cpp
 issm_CXXFLAGS= $(CXXFLAGS) $(CXXOPTFLAGS)
 
Index: /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.cpp	(revision 26000)
+++ /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.cpp	(revision 26000)
@@ -0,0 +1,420 @@
+#include "./SamplingAnalysis.h" //
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include <random>
+
+#define FINITEELEMENT P1Enum
+#define 	NUMVERTICES   3
+
+/* Reimplementation of IsFaceOnBoundary with two edges on boundary supported */
+bool IsFaceOnBoundary(Element*  element)
+{
+
+	IssmDouble values[NUMVERTICES];
+	IssmDouble sum;
+
+	/*Retrieve all inputs and parameters*/
+	element->GetInputListOnVertices(&values[0],MeshVertexonboundaryEnum);
+	sum = values[0]+values[1]+values[2];
+
+	_assert_(sum==0. || sum==1. || sum==2. || sum==3.);
+
+	if(sum>1.){
+	 	return true;
+	}
+	else{
+	 return false;
+	}
+}
+
+/*Model processing*/
+void SamplingAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+	/*Default, do nothing*/
+	return;
+}/*}}}*/
+void SamplingAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
+	/*Default, do nothing*/
+	return;
+}/*}}}*/
+void SamplingAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
+
+	int finiteelement;
+	finiteelement = FINITEELEMENT;
+
+	/*Check in 2d*/
+	if(iomodel->domaintype!=Domain2DhorizontalEnum) _error_("Only 2D horizontal domain is implemented.");
+
+	::CreateNodes(nodes,iomodel,SamplingAnalysisEnum,finiteelement);
+
+}/*}}}*/
+int  SamplingAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
+	return 1;
+}/*}}}*/
+void SamplingAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
+
+	int    finiteelement;
+
+	/*Finite element type*/
+	finiteelement = FINITEELEMENT;
+
+	/*Update elements: */
+	int counter=0;
+	for(int i=0;i<iomodel->numberofelements;i++){
+		if(iomodel->my_elements[i]){
+			Element* element=(Element*)elements->GetObjectByOffset(counter);
+			element->Update(inputs,i,iomodel,analysis_counter,analysis_type,finiteelement);
+			counter++;
+		}
+	}
+
+	/*Create inputs: */
+	iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
+	iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonboundary",MeshVertexonboundaryEnum);
+
+  iomodel->FetchDataToInput(inputs,elements,"md.sampling.kappa",SamplingKappaEnum);
+	iomodel->FetchDataToInput(inputs,elements,"md.sampling.beta",SamplingBetaEnum,0.);
+	iomodel->FetchDataToInput(inputs,elements,"md.initialization.sample",SampleEnum,0.);
+
+}/*}}}*/
+void SamplingAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
+
+	int     numoutputs;
+	char**  requestedoutputs = NULL;
+
+	parameters->AddObject(iomodel->CopyConstantObject("md.sampling.tau",SamplingTauEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.sampling.alpha",SamplingAlphaEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.sampling.robin",SamplingRobinEnum));
+	if(solution_enum==TransientSolutionEnum) parameters->AddObject(iomodel->CopyConstantObject("md.sampling.phi",SamplingPhiEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.sampling.seed",SamplingSeedEnum));
+
+	iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.sampling.requested_outputs");
+	parameters->AddObject(new IntParam(SamplingNumRequestedOutputsEnum,numoutputs));
+	if(numoutputs)parameters->AddObject(new StringArrayParam(SamplingRequestedOutputsEnum,requestedoutputs,numoutputs));
+	iomodel->DeleteData(&requestedoutputs,numoutputs,"md.sampling.requested_outputs");
+
+}/*}}}*/
+
+/*Finite Element Analysis*/
+void           SamplingAnalysis::Core(FemModel* femmodel){/*{{{*/
+	_error_("not implemented");
+}/*}}}*/
+ElementVector* SamplingAnalysis::CreateDVector(Element* element){/*{{{*/
+	/*Default, return NULL*/
+	return NULL;
+}/*}}}*/
+ElementMatrix* SamplingAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
+_error_("not implemented");
+}/*}}}*/
+ElementMatrix* SamplingAnalysis::CreateKMatrix(Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	//if(!element->IsIceInElement()) return NULL; // Ice in element not required for sampling capability
+
+	/*Intermediaries*/
+	bool robin;
+	int      domaintype;
+
+	/*compute all stiffness matrices for this element*/
+
+	ElementMatrix* Ke = NULL;
+	element->FindParam(&robin,SamplingRobinEnum);
+
+  if(!robin) Ke=CreateKMatrixModifiedHelmholtz(element);
+  else{
+		ElementMatrix* Ke1=CreateKMatrixModifiedHelmholtz(element);
+		ElementMatrix* Ke2=CreateKMatrixRobinBC(element);
+		Ke =new ElementMatrix(Ke1,Ke2);
+
+		delete Ke1;
+		delete Ke2;
+  }
+
+	/*clean-up and return*/
+	return Ke;
+
+}/*}}}*/
+ElementMatrix* SamplingAnalysis::CreateKMatrixModifiedHelmholtz(Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	//if(!element->IsIceInElement()) return NULL;
+
+	/*Intermediaries*/
+	int dim;
+	IssmDouble  D,Jdet;
+	IssmDouble* xyz_list = NULL;
+	IssmDouble kappa;
+
+	/*Get problem dimension*/
+	dim=2;	// Only 2D horizontal domain is implemented so far
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementMatrix* Ke     = element->NewElementMatrix();
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input2*	kappa_input=element->GetInput2(SamplingKappaEnum); _assert_(kappa_input);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+		kappa_input->GetInputValue(&kappa,gauss); _assert_(kappa>0);
+
+		D=gauss->weight*Jdet;
+
+		/* Laplacian operator */
+		for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j] += D*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]);
+			}
+ 		}
+
+		/* Identity identity */
+		for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j] += D*kappa*kappa*(basis[j]*basis[i]);
+			}
+ 		}
+
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
+	delete gauss;
+	return Ke;
+
+}/*}}}*/
+ElementMatrix* SamplingAnalysis::CreateKMatrixRobinBC(Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	//if(!element->IsIceInElement()) return NULL;
+
+	/*If no boundary, return NULL*/
+	if(!IsFaceOnBoundary(element)) return NULL;
+
+	/*Intermediaries*/
+	IssmDouble  D,Jdet;
+	IssmDouble* xyz_list = NULL;
+	IssmDouble* xyz_list_boundary = NULL;
+	IssmDouble beta;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementMatrix* Ke     = element->NewElementMatrix();
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->GetLevelCoordinates(&xyz_list_boundary,xyz_list,MeshVertexonboundaryEnum,1.);
+	Input2*	beta_input=element->GetInput2(SamplingBetaEnum); _assert_(beta_input);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(xyz_list,xyz_list_boundary,3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminantSurface(&Jdet,xyz_list_boundary,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		beta_input->GetInputValue(&beta,gauss); _assert_(beta>=0);
+
+		D=gauss->weight*Jdet;
+
+		for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j] += D*beta*(basis[j]*basis[i]);
+			}
+ 		}
+
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(xyz_list_boundary);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return Ke;
+
+}/*}}}*/
+ElementVector* SamplingAnalysis::CreatePVector(Element* element){/*{{{*/
+	//_error_("not supported");
+	return NULL;
+}/*}}}*/
+void           SamplingAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+	element->GetSolutionFromInputsOneDof(solution,SampleEnum);
+}/*}}}*/
+void           SamplingAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
+	_error_("not supported");
+}/*}}}*/
+void           SamplingAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+
+    /*Fetch number of nodes and dof for this finite element*/
+    int numnodes = element->GetNumberOfNodes();
+
+		/*Fetch dof list and allocate solution vector*/
+		int* doflist=NULL;
+    element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
+    IssmDouble* newsample = xNew<IssmDouble>(numnodes);
+
+    /*Use the dof list to index into the solution vector: */
+    for(int i=0;i<numnodes;i++){
+			newsample[i]=solution[doflist[i]];
+			/*Check solution*/
+			if(xIsNan<IssmDouble>(newsample[i])) _error_("NaN found in solution vector");
+			if(xIsInf<IssmDouble>(newsample[i])) _error_("Inf found in solution vector");
+		}
+
+    /*Add sample inputs to the tria element: */
+    element->AddInput2(SampleEnum,newsample,element->GetElementType());
+
+    /*Free ressources:*/
+    xDelete<IssmDouble>(newsample);
+    xDelete<int>(doflist);
+ }/*}}}*/
+void           SamplingAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+	/*Default, do nothing*/
+  return;
+}/*}}}*/
+
+ElementMatrix* SamplingAnalysis::CreateMassMatrix(Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	//if(!element->IsIceInElement()) return NULL;
+
+	/*Intermediaries*/
+	IssmDouble  D,Jdet;
+	IssmDouble* xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementMatrix* Me     = element->NewElementMatrix();
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		D=gauss->weight*Jdet;
+		TripleMultiply(basis,1,numnodes,1,
+					&D,1,1,0,
+					basis,1,numnodes,0,
+					&Me->values[0],1);
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return Me;
+}/*}}}*/
+void           SamplingAnalysis::LumpedKMatrix(Vector<IssmDouble>** pKlff,FemModel* femmodel){/*{{{*/
+
+	/*Initialize Lumped mass matrix (actually we just save its diagonal)*/
+	int fsize      = femmodel->nodes->NumberOfDofs(FsetEnum);
+	int flocalsize = femmodel->nodes->NumberOfDofsLocal(FsetEnum);
+	Vector<IssmDouble>* Klff = new Vector<IssmDouble>(flocalsize,fsize);
+
+	/*Create and assemble matrix*/
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element*       element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		ElementMatrix* KLe     = this->CreateKMatrix(element);
+		if(KLe){
+			KLe->Lump();
+			KLe->AddDiagonalToGlobal(Klff);
+		}
+		delete KLe;
+	}
+	Klff->Assemble();
+
+	/*Assign output pointer*/
+	*pKlff=Klff;
+}/*}}}*/
+void           SamplingAnalysis::KMatrix(Matrix<IssmDouble>** pKff,FemModel* femmodel){/*{{{*/
+
+	/*Initialize K matrix*/
+	Matrix<IssmDouble> *Kff = NULL;
+	AllocateSystemMatricesx(&Kff,NULL,NULL,NULL,femmodel);
+
+	/*Create and assemble matrix*/
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element*       element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		ElementMatrix* Ke     = this->CreateKMatrix(element);
+		if(Ke){
+			Ke->AddToGlobal(Kff);
+		}
+		delete Ke;
+	}
+	Kff->Assemble();
+
+	/*Assign output pointer*/
+	*pKff=Kff;
+
+}/*}}}*/
+void           SamplingAnalysis::LumpedMassMatrix(Vector<IssmDouble>** pMlff,FemModel* femmodel){/*{{{*/
+
+	/*Initialize Lumped mass matrix (actually we just save its diagonal)*/
+	int fsize      = femmodel->nodes->NumberOfDofs(FsetEnum);
+	int flocalsize = femmodel->nodes->NumberOfDofsLocal(FsetEnum);
+	Vector<IssmDouble>* Mlff = new Vector<IssmDouble>(flocalsize,fsize);
+
+	/*Create and assemble matrix*/
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element*       element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		ElementMatrix* MLe     = this->CreateMassMatrix(element);
+		if(MLe){
+			MLe->Lump();
+			MLe->AddDiagonalToGlobal(Mlff);
+		}
+		delete MLe;
+	}
+	Mlff->Assemble();
+
+	/*Assign output pointer*/
+	*pMlff=Mlff;
+}/*}}}*/
+void           SamplingAnalysis::MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel){/*{{{*/
+
+	/*Initialize Mass matrix*/
+	Matrix<IssmDouble> *Mff = NULL;
+	AllocateSystemMatricesx(&Mff,NULL,NULL,NULL,femmodel);
+
+	/*Create and assemble matrix*/
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element*       element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		ElementMatrix* MLe     = this->CreateMassMatrix(element);
+		if(MLe){
+			MLe->AddToGlobal(Mff);
+		}
+		delete MLe;
+	}
+	Mff->Assemble();
+
+	/*Assign output pointer*/
+	*pMff=Mff;
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.h	(revision 26000)
+++ /issm/trunk-jpl/src/c/analyses/SamplingAnalysis.h	(revision 26000)
@@ -0,0 +1,42 @@
+/*! \file SamplingAnalysis.h
+ *  \brief: header file for generic external result object
+ */
+
+#ifndef _SamplingAnalysis_
+#define _SamplingAnalysis_
+
+/*Headers*/
+#include "./Analysis.h"
+
+class SamplingAnalysis: public Analysis{
+
+	public:
+		/*Model processing*/
+		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
+		void CreateLoads(Loads* loads, IoModel* iomodel);
+		void CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr=false);
+		int  DofsPerNode(int** doflist,int domaintype,int approximation);
+		void UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type);
+		void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
+
+		/*Finite element Analysis*/
+		void           Core(FemModel* femmodel);
+		ElementVector* CreateDVector(Element* element);
+		ElementMatrix* CreateJacobianMatrix(Element* element);
+		ElementMatrix* CreateKMatrix(Element* element);
+		ElementMatrix* CreateKMatrixModifiedHelmholtz(Element* element);
+		ElementMatrix* CreateKMatrixRobinBC(Element* element);
+		ElementMatrix* CreateMassMatrix(Element* element);
+    ElementVector* CreatePVector(Element* element);
+		void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
+		void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
+		void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
+		void           UpdateConstraints(FemModel* femmodel);
+
+		void           LumpedKMatrix(Vector<IssmDouble>** pKLff,FemModel* femmodel);
+		void           KMatrix(Matrix<IssmDouble>** pKff,FemModel* femmodel);
+		void           LumpedMassMatrix(Vector<IssmDouble>** pMLff,FemModel* femmodel);
+		void           MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel);
+
+};
+#endif
Index: /issm/trunk-jpl/src/c/analyses/analyses.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/analyses.h	(revision 25999)
+++ /issm/trunk-jpl/src/c/analyses/analyses.h	(revision 26000)
@@ -1,4 +1,4 @@
 /*
- * analyses.h: 
+ * analyses.h:
  */
 
@@ -35,4 +35,5 @@
 #include "./LevelsetAnalysis.h"
 #include "./MasstransportAnalysis.h"
+#include "./SamplingAnalysis.h"
 #include "./SmbAnalysis.h"
 #include "./SealevelchangeAnalysis.h"