Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 17830)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 17831)
@@ -430,4 +430,5 @@
 					./solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp\
 					./cores/balancethickness_core.cpp \
+					./cores/balancethickness2_core.cpp \
 					./cores/balancevelocity_core.cpp \
 					./cores/dummy_core.cpp\
@@ -521,4 +522,7 @@
 if BALANCETHICKNESS
 issm_sources += ./analyses/BalancethicknessAnalysis.cpp
+endif
+if BALANCETHICKNESS2
+issm_sources += ./analyses/Balancethickness2Analysis.cpp
 endif
 if BALANCETHICKNESSSOFT
@@ -866,5 +870,5 @@
 
 libISSMCore_la_SOURCES  = $(issm_sources)
-libISSMCore_la_LIBADD = $(PETSCLIB) $(TAOLIB) $(PLAPACKLIB) $(MUMPSLIB) $(SUPERLULIB) $(SPOOLESLIB) $(SCALAPACKLIB) $(BLACSLIB) $(HYPRELIB) $(SPAILIB) $(PROMETHEUSLIB) $(PASTIXLIB) $(MLLIB) $(DAKOTALIB) $(METISLIB) $(CHACOLIB) $(SCOTCHLIB) $(BLASLAPACKLIB) $(MKLLIB) $(MPILIB) $(MATHLIB) $(FORTRANLIB) $(GRAPHICSLIB) $(MULTITHREADINGLIB) $(OSLIBS) $(GSLLIB)   $(ADOLCLIB) $(AMPILIB) $(METEOIOLIB) $(SNOWPACKLIB) 
+libISSMCore_la_LIBADD = $(PETSCLIB) $(TAOLIB) $(PLAPACKLIB) $(MUMPSLIB) $(SUPERLULIB) $(SPOOLESLIB) $(SCALAPACKLIB) $(BLACSLIB) $(HYPRELIB) $(SPAILIB) $(PROMETHEUSLIB) $(PASTIXLIB) $(MLLIB) $(DAKOTALIB) $(METISLIB) $(CHACOLIB) $(SCOTCHLIB) $(BLASLAPACKLIB) $(MKLLIB) $(MPILIB) $(MATHLIB) $(GRAPHICSLIB) $(MULTITHREADINGLIB) $(OSLIBS) $(GSLLIB)   $(ADOLCLIB) $(AMPILIB) $(METEOIOLIB) $(SNOWPACKLIB) 
 libISSMCore_la_FFLAGS = $(AM_FFLAGS)
 endif
@@ -913,8 +917,8 @@
 
 #External packages
-LDADD += $(PETSCLIB) $(TAOLIB) $(PLAPACKLIB) $(MUMPSLIB) $(SUPERLULIB) $(SPOOLESLIB) $(SCALAPACKLIB) $(BLACSLIB) $(HYPRELIB) $(SPAILIB) $(PROMETHEUSLIB) $(PASTIXLIB) $(MLLIB) $(DAKOTALIB) $(METISLIB) $(CHACOLIB) $(SCOTCHLIB) $(BLASLAPACKLIB) $(MKLLIB) $(MPILIB)  $(MATHLIB) $(FORTRANLIB) $(GRAPHICSLIB) $(MULTITHREADINGLIB) $(OSLIBS) $(GSLLIB) $(AMPILIB) $(ADOLCLIB) $(MPILIB) $(METEOIOLIB) $(SNOWPACKLIB)
+LDADD += $(PETSCLIB) $(TAOLIB) $(PLAPACKLIB) $(MUMPSLIB) $(SUPERLULIB) $(SPOOLESLIB) $(SCALAPACKLIB) $(BLACSLIB) $(HYPRELIB) $(SPAILIB) $(PROMETHEUSLIB) $(PASTIXLIB) $(MLLIB) $(DAKOTALIB) $(METISLIB) $(CHACOLIB) $(SCOTCHLIB) $(BLASLAPACKLIB) $(MKLLIB) $(MPILIB)  $(MATHLIB) $(GRAPHICSLIB) $(MULTITHREADINGLIB) $(OSLIBS) $(GSLLIB) $(AMPILIB) $(ADOLCLIB) $(MPILIB) $(METEOIOLIB) $(SNOWPACKLIB)
 
 if FORTRAN
-LDADD += $(FLIBS)
+LDADD += $(FLIBS) $(FORTRANLIB)
 endif
 
Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 17831)
+++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 17831)
@@ -0,0 +1,304 @@
+#include "./Balancethickness2Analysis.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+
+/*Model processing*/
+int  Balancethickness2Analysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
+	return 1;
+}/*}}}*/
+void Balancethickness2Analysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
+}/*}}}*/
+void Balancethickness2Analysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
+
+	/*Finite element type*/
+	int finiteelement = P1Enum;
+
+	/*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(i,iomodel,analysis_counter,analysis_type,finiteelement);
+			counter++;
+		}
+	}
+
+	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
+	iomodel->FetchDataToInput(elements,VxEnum);
+	iomodel->FetchDataToInput(elements,VyEnum);
+	iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
+	iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
+	iomodel->FetchDataToInput(elements,BalancethicknessThickeningRateEnum);
+	iomodel->FetchDataToInput(elements,MeshVertexonboundaryEnum);
+	iomodel->FetchDataToInput(elements,ThicknessEnum);
+}/*}}}*/
+void Balancethickness2Analysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
+
+	int finiteelement = P1Enum;
+	::CreateNodes(nodes,iomodel,Balancethickness2AnalysisEnum,finiteelement);
+}/*}}}*/
+void Balancethickness2Analysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+
+	/*Do not add constraints in DG*/
+	IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,Balancethickness2AnalysisEnum,P1Enum);
+
+}/*}}}*/
+void Balancethickness2Analysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
+
+}/*}}}*/
+
+/*Finite Element Analysis*/
+void           Balancethickness2Analysis::Core(FemModel* femmodel){/*{{{*/
+	_error_("not implemented yet");
+}/*}}}*/
+ElementVector* Balancethickness2Analysis::CreateDVector(Element* element){/*{{{*/
+	/*Default, return NULL*/
+	return NULL;
+}/*}}}*/
+ElementMatrix* Balancethickness2Analysis::CreateJacobianMatrix(Element* element){/*{{{*/
+_error_("Not implemented");
+}/*}}}*/
+ElementMatrix* Balancethickness2Analysis::CreateKMatrix(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	IssmDouble Jdet,D_scalar;
+	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* Ke  = element->NewElementMatrix();
+	IssmDouble*    B   = xNew<IssmDouble>(2*numnodes);
+	IssmDouble     D[2][2];
+
+	/*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);
+
+		IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
+			}
+		}
+		//GetB(B,element,xyz_list,gauss);
+
+		//D_scalar=gauss->weight*Jdet;
+
+		//D[0][0]=D_scalar;
+		//D[0][1]=0.;
+		//D[1][0]=0.;
+		//D[1][1]=D_scalar;
+		//TripleMultiply(B,2,numnodes,1,
+		//			&D[0][0],2,2,0,
+		//			B,2,numnodes,0,
+		//			&Ke->values[0],1);
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(B);
+	delete gauss;
+	return Ke;
+}/*}}}*/
+ElementVector* Balancethickness2Analysis::CreatePVector(Element* element){/*{{{*/
+
+	/*compute all load vectors for this element*/
+	ElementVector* pe1=CreatePVectorVolume(element);
+	ElementVector* pe2=CreatePVectorBoundary(element);
+	ElementVector* pe =new ElementVector(pe1,pe2);
+
+	/*clean-up and return*/
+	delete pe1;
+	delete pe2;
+	return pe;
+}/*}}}*/
+ElementVector* Balancethickness2Analysis::CreatePVectorVolume(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	IssmDouble  dhdt,mb,ms,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*/
+	ElementVector* pe    = element->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input* mb_input   = element->GetInput(BasalforcingsMeltingRateEnum);       _assert_(mb_input);
+	Input* ms_input   = element->GetInput(SurfaceforcingsMassBalanceEnum);     _assert_(ms_input);
+	Input* dhdt_input = element->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
+
+	/*Initialize mb_correction to 0, do not forget!:*/
+	/* 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);
+
+		ms_input->GetInputValue(&ms,gauss);
+		mb_input->GetInputValue(&mb,gauss);
+		dhdt_input->GetInputValue(&dhdt,gauss);
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms-mb-dhdt)*basis[i];
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return pe;
+}/*}}}*/
+ElementVector* Balancethickness2Analysis::CreatePVectorBoundary(Element* element){/*{{{*/
+
+	/*If no front, return NULL*/
+	if(!element->IsFaceOnBoundary()) return NULL;
+
+	/*Intermediaries*/
+	IssmDouble  Jdet,thickness,vx,vy;
+	IssmDouble *xyz_list = NULL;
+	IssmDouble *xyz_list_front = NULL;
+	IssmDouble  normal[2];
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = element->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
+	Input* vx_input        = element->GetInput(VxEnum);        _assert_(vx_input);
+	Input* vy_input        = element->GetInput(VyEnum);        _assert_(vy_input);
+
+	element->GetVerticesCoordinates(&xyz_list);
+	element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
+	element->NormalSection(&normal[0],xyz_list_front);
+
+	/*Start looping on Gaussian points*/
+	Gauss* gauss=element->NewGauss(xyz_list,xyz_list_front,3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+		thickness_input->GetInputValue(&thickness,gauss);
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		for(int i=0;i<numnodes;i++) pe->values[i] += Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1])*basis[i];
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(xyz_list_front);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return pe;
+}/*}}}*/
+void Balancethickness2Analysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	 * For node i, Bi can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi=[ dN/dx ]
+	 *          [ dN/dy ]
+	 * where N is the finiteelement function for node i.
+	 *
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+	/*Build B': */
+	for(int i=0;i<numnodes;i++){
+		B[numnodes*0+i] = dbasis[0*numnodes+i];
+		B[numnodes*1+i] = dbasis[1*numnodes+i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
+
+}/*}}}*/
+void Balancethickness2Analysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+	   _error_("not implemented yet");
+}/*}}}*/
+void Balancethickness2Analysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+
+	/*Intermediaries */
+	IssmDouble  vx,vy,vel,h,normdphi,dphi[2];
+	IssmDouble* xyz_list = NULL;
+	int       * doflist  = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes    = element->GetNumberOfNodes();
+	int numvertices = element->GetNumberOfVertices();
+
+	/*Fetch dof list and allocate solution vector*/
+	element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	IssmDouble* thickness = xNew<IssmDouble>(numvertices);
+	IssmDouble* values    = xNew<IssmDouble>(numnodes);
+
+	/*Use the dof list to index into the solution vector: */
+	for(int i=0;i<numnodes;i++){
+		values[i]=solution[doflist[i]];
+
+		/*Check solution*/
+		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
+	}
+
+	element->AddInput(EnthalpyEnum,values,element->GetElementType());//FIXME
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input* potential_input = element->GetInput(EnthalpyEnum); _assert_(potential_input);
+	Input* vx_input        = element->GetInput(VxEnum);       _assert_(vx_input);
+	Input* vy_input        = element->GetInput(VyEnum);       _assert_(vy_input);
+
+	Gauss* gauss=element->NewGauss();
+	for (int iv=0;iv<numvertices;iv++){
+		gauss->GaussVertex(iv);
+
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		potential_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
+
+		vel = sqrt(vx*vx + vy*vy) + 1.e-10;
+		normdphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]);
+
+		thickness[iv] = normdphi/vel;
+	}
+	element->AddInput(ThicknessEnum,thickness,P1Enum);
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<int>(doflist);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(values);
+	xDelete<IssmDouble>(thickness);
+}/*}}}*/
+void Balancethickness2Analysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+	/*Default, do nothing*/
+	return;
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h	(revision 17831)
+++ /issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h	(revision 17831)
@@ -0,0 +1,35 @@
+/*! \file Balancethickness2Analysis.h 
+ *  \brief: header file for generic external result object
+ */
+
+#ifndef _Balancethickness2Analysis_
+#define _Balancethickness2Analysis_
+
+/*Headers*/
+#include "./Analysis.h"
+
+class Balancethickness2Analysis: public Analysis{
+
+	public:
+		/*Model processing*/
+		int  DofsPerNode(int** doflist,int domaintype,int approximation);
+		void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
+		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
+		void CreateNodes(Nodes* nodes,IoModel* iomodel);
+		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
+		void CreateLoads(Loads* loads, IoModel* iomodel);
+
+		/*Finite element Analysis*/
+		void           Core(FemModel* femmodel);
+		ElementVector* CreateDVector(Element* element);
+		ElementMatrix* CreateJacobianMatrix(Element* element);
+		ElementMatrix* CreateKMatrix(Element* element);
+		ElementVector* CreatePVector(Element* element);
+		ElementVector* CreatePVectorVolume(Element* element);
+		ElementVector* CreatePVectorBoundary(Element* element);
+		void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
+		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
+		void UpdateConstraints(FemModel* femmodel);
+};
+#endif
Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp	(revision 17830)
+++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp	(revision 17831)
@@ -22,4 +22,7 @@
 		#ifdef _HAVE_BALANCETHICKNESS_
 		case BalancethicknessAnalysisEnum : return new BalancethicknessAnalysis();
+		#endif
+		#ifdef _HAVE_BALANCETHICKNESS2_
+		case Balancethickness2AnalysisEnum : return new Balancethickness2Analysis();
 		#endif
 		#ifdef _HAVE_BALANCETHICKNESSSOFT_
Index: /issm/trunk-jpl/src/c/analyses/analyses.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/analyses.h	(revision 17830)
+++ /issm/trunk-jpl/src/c/analyses/analyses.h	(revision 17831)
@@ -11,4 +11,5 @@
 #include "./AdjointHorizAnalysis.h"
 #include "./BalancethicknessAnalysis.h"
+#include "./Balancethickness2Analysis.h"
 #include "./BalancethicknessSoftAnalysis.h"
 #include "./BalancevelocityAnalysis.h"
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17830)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17831)
@@ -243,4 +243,5 @@
 		virtual bool   IsZeroLevelset(int levelset_enum)=0;
 		virtual bool   IsIcefront(void)=0;
+		virtual bool   IsFaceOnBoundary(void)=0;
 		virtual void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum)=0;
 		virtual void   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17830)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17831)
@@ -95,4 +95,5 @@
 		bool   IsZeroLevelset(int levelset_enum);
 		bool   IsIcefront(void);
+		bool   IsFaceOnBoundary(void){_error_("not implemented yet");};
 		void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
 		void   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17830)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17831)
@@ -126,4 +126,5 @@
 		bool        IsZeroLevelset(int levelset_enum){_error_("not implemented");};
 		bool		   IsIcefront(void);
+		bool   IsFaceOnBoundary(void){_error_("not implemented yet");};
 		void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");};
 		void		   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 17830)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 17831)
@@ -131,4 +131,5 @@
 		bool        IsZeroLevelset(int levelset_enum){_error_("not implemented");};
 		bool		   IsIcefront(void);
+		bool   IsFaceOnBoundary(void){_error_("not implemented yet");};
 		void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
 		void		   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17830)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17831)
@@ -2257,4 +2257,24 @@
 	return isicefront;
 }/*}}}*/
+bool Tria::IsFaceOnBoundary(void){/*{{{*/
+
+	IssmDouble values[NUMVERTICES];
+	IssmDouble sum;
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&values[0],MeshVertexonboundaryEnum);
+	sum = values[0]+values[1]+values[2];
+
+	_assert_(sum==0. || sum==1. || sum==2.);
+
+	if(sum==3.)  _error_("Two edges on boundary not supported yet...");
+
+	if(sum>1.){
+		return true;
+	}
+	else{
+		return false;
+	}
+}/*}}}*/
 /*FUNCTION Tria::AverageOntoPartition {{{*/
 void  Tria::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17830)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17831)
@@ -112,4 +112,5 @@
 		bool        IsZeroLevelset(int levelset_enum);
 		bool		IsIcefront(void);
+		bool		IsFaceOnBoundary(void);
 
 		void       AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp	(revision 17830)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp	(revision 17831)
@@ -31,5 +31,4 @@
 	this->helement=NULL;
 	this->element=NULL;
-	this->isdamaged=NULL;
 	return;
 }
Index: /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp	(revision 17830)
+++ /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp	(revision 17831)
@@ -78,4 +78,10 @@
 			break;
 
+		case Balancethickness2SolutionEnum:
+			numanalyses=1;
+			analyses=xNew<int>(numanalyses);
+			analyses[0]=Balancethickness2AnalysisEnum;
+			break;
+
 		case BalancethicknessSoftSolutionEnum:
 			numanalyses=1;
Index: /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp	(revision 17830)
+++ /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp	(revision 17831)
@@ -34,4 +34,7 @@
 		case BalancethicknessSolutionEnum:
 			solutioncore=&balancethickness_core;
+			break;
+		case Balancethickness2SolutionEnum:
+			solutioncore=&balancethickness2_core;
 			break;
 		case BalancethicknessSoftSolutionEnum:
Index: /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp	(revision 17831)
+++ /issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp	(revision 17831)
@@ -0,0 +1,32 @@
+/*!\file: balancethickness_core.cpp
+ * \brief: core of the balancethickness solution 
+ */ 
+
+#include "./cores.h"
+#include "../toolkits/toolkits.h"
+#include "../classes/classes.h"
+#include "../shared/shared.h"
+#include "../modules/modules.h"
+#include "../solutionsequences/solutionsequences.h"
+
+void balancethickness2_core(FemModel* femmodel){
+
+	/*parameters: */
+	bool save_results;
+
+	/*activate formulation: */
+	femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
+
+	/*recover parameters: */
+	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+
+	if(VerboseSolution()) _printf0_("call computational core:\n");
+	solutionsequence_linear(femmodel);
+
+	if(save_results){
+		if(VerboseSolution()) _printf0_("   saving results\n");
+		int outputs[2] = {ThicknessEnum,EnthalpyEnum};
+		femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
+	}
+
+}
Index: /issm/trunk-jpl/src/c/cores/cores.h
===================================================================
--- /issm/trunk-jpl/src/c/cores/cores.h	(revision 17830)
+++ /issm/trunk-jpl/src/c/cores/cores.h	(revision 17831)
@@ -34,4 +34,5 @@
 void extrudefromtop_core(FemModel* femmodel);
 void balancethickness_core(FemModel* femmodel);
+void balancethickness2_core(FemModel* femmodel);
 void balancevelocity_core(FemModel* femmodel);
 void slopecompute_core(FemModel* femmodel);
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 17830)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 17831)
@@ -212,4 +212,5 @@
 	MeshVertexonbaseEnum,
 	MeshVertexonsurfaceEnum,
+	MeshVertexonboundaryEnum,
 	MeshXEnum,
 	MeshYEnum,
@@ -324,4 +325,6 @@
 	BalancethicknessAnalysisEnum,
 	BalancethicknessSolutionEnum,
+	Balancethickness2AnalysisEnum,
+	Balancethickness2SolutionEnum,
 	BalancethicknessSoftAnalysisEnum,
 	BalancethicknessSoftSolutionEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 17830)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 17831)
@@ -220,4 +220,5 @@
 		case MeshVertexonbaseEnum : return "MeshVertexonbase";
 		case MeshVertexonsurfaceEnum : return "MeshVertexonsurface";
+		case MeshVertexonboundaryEnum : return "MeshVertexonboundary";
 		case MeshXEnum : return "MeshX";
 		case MeshYEnum : return "MeshY";
@@ -325,4 +326,6 @@
 		case BalancethicknessAnalysisEnum : return "BalancethicknessAnalysis";
 		case BalancethicknessSolutionEnum : return "BalancethicknessSolution";
+		case Balancethickness2AnalysisEnum : return "Balancethickness2Analysis";
+		case Balancethickness2SolutionEnum : return "Balancethickness2Solution";
 		case BalancethicknessSoftAnalysisEnum : return "BalancethicknessSoftAnalysis";
 		case BalancethicknessSoftSolutionEnum : return "BalancethicknessSoftSolution";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 17830)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 17831)
@@ -223,4 +223,5 @@
 	      else if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum;
 	      else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
+	      else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum;
 	      else if (strcmp(name,"MeshX")==0) return MeshXEnum;
 	      else if (strcmp(name,"MeshY")==0) return MeshYEnum;
@@ -259,9 +260,9 @@
 	      else if (strcmp(name,"SettingsLowmem")==0) return SettingsLowmemEnum;
 	      else if (strcmp(name,"SettingsOutputFrequency")==0) return SettingsOutputFrequencyEnum;
-	      else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum;
          else stage=3;
    }
    if(stage==3){
-	      if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum;
+	      if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum;
+	      else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum;
 	      else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum;
 	      else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum;
@@ -331,4 +332,6 @@
 	      else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
 	      else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
+	      else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
+	      else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
 	      else if (strcmp(name,"BalancethicknessSoftAnalysis")==0) return BalancethicknessSoftAnalysisEnum;
 	      else if (strcmp(name,"BalancethicknessSoftSolution")==0) return BalancethicknessSoftSolutionEnum;
@@ -380,11 +383,11 @@
 	      else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
 	      else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
-	      else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum;
-	      else if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum;
-	      else if (strcmp(name,"HOFSApproximation")==0) return HOFSApproximationEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
+	      if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum;
+	      else if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum;
+	      else if (strcmp(name,"HOFSApproximation")==0) return HOFSApproximationEnum;
+	      else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
 	      else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum;
 	      else if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
@@ -503,11 +506,11 @@
 	      else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum;
 	      else if (strcmp(name,"Temperature")==0) return TemperatureEnum;
-	      else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
-	      else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
-	      else if (strcmp(name,"Vel")==0) return VelEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"Velocity")==0) return VelocityEnum;
+	      if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
+	      else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
+	      else if (strcmp(name,"Vel")==0) return VelEnum;
+	      else if (strcmp(name,"Velocity")==0) return VelocityEnum;
 	      else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
 	      else if (strcmp(name,"Vx")==0) return VxEnum;
@@ -626,11 +629,11 @@
 	      else if (strcmp(name,"AugmentedLagrangianTheta")==0) return AugmentedLagrangianThetaEnum;
 	      else if (strcmp(name,"None")==0) return NoneEnum;
-	      else if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum;
-	      else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum;
-	      else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"SubelementMigration2")==0) return SubelementMigration2Enum;
+	      if (strcmp(name,"AgressiveMigration")==0) return AgressiveMigrationEnum;
+	      else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum;
+	      else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
+	      else if (strcmp(name,"SubelementMigration2")==0) return SubelementMigration2Enum;
 	      else if (strcmp(name,"Contact")==0) return ContactEnum;
 	      else if (strcmp(name,"MaskGroundediceLevelset")==0) return MaskGroundediceLevelsetEnum;
Index: /issm/trunk-jpl/src/c/shared/Enum/Synchronize.sh
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Synchronize.sh	(revision 17830)
+++ /issm/trunk-jpl/src/c/shared/Enum/Synchronize.sh	(revision 17831)
@@ -35,5 +35,5 @@
 END
 #core
-cat temp | grep [a-zA-Z]Analysis | \
+cat temp | grep [a-zA-Z0-9]Analysis | \
 	grep -v DefaultAnalysis | grep -v FlaimAnalysis | grep -v SurfaceSlopeAnalysis | grep -v BedSlopeAnalysis | \
 	awk '{print "\t\t#ifdef _HAVE_"toupper(substr($1,1,length($1)-12))"_\n\t\t" "case " $1" : return new " substr($1,1,length($1)-4) "();\n\t\t#endif"}' \
@@ -63,5 +63,5 @@
 END
 #core
-cat temp | grep [a-zA-Z]Analysis | \
+cat temp | grep [a-zA-Z0-9]Analysis | \
 	grep -v DefaultAnalysis | grep -v FlaimAnalysis | grep -v SurfaceSlopeAnalysis | grep -v BedSlopeAnalysis | \
 	sed -e "s/AnalysisEnum//g" | \
Index: /issm/trunk-jpl/src/m/classes/mesh2d.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/mesh2d.m	(revision 17830)
+++ /issm/trunk-jpl/src/m/classes/mesh2d.m	(revision 17831)
@@ -190,4 +190,5 @@
 			WriteData(fid,'object',obj,'class','mesh','fieldname','numberofvertices','format','Integer');
 			WriteData(fid,'object',obj,'class','mesh','fieldname','average_vertex_connectivity','format','Integer');
+			WriteData(fid,'object',obj,'class','mesh','fieldname','vertexonboundary','format','DoubleMat','mattype',1);
 		end % }}}
 		function t = domaintype(obj) % {{{
Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m
===================================================================
--- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m	(revision 17830)
+++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m	(revision 17831)
@@ -61,4 +61,7 @@
 		analyses=[BalancethicknessAnalysisEnum()];
 
+	case Balancethickness2SolutionEnum(),
+		analyses=[Balancethickness2AnalysisEnum()];
+
 	case BalancethicknessSoftSolutionEnum(),
 		analyses=[BalancethicknessAnalysisEnum()];
Index: /issm/trunk-jpl/src/m/contrib/dassflow/importgmsh.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/dassflow/importgmsh.m	(revision 17830)
+++ /issm/trunk-jpl/src/m/contrib/dassflow/importgmsh.m	(revision 17831)
@@ -97,4 +97,6 @@
 				error('not supported');
 			end
+		case 15, %point
+			continue;
 		otherwise,
 			error(['Type ' num2str(ty) ' not supported']);
Index: /issm/trunk-jpl/src/m/enum/Balancethickness2AnalysisEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/Balancethickness2AnalysisEnum.m	(revision 17831)
+++ /issm/trunk-jpl/src/m/enum/Balancethickness2AnalysisEnum.m	(revision 17831)
@@ -0,0 +1,11 @@
+function macro=Balancethickness2AnalysisEnum()
+%BALANCETHICKNESS2ANALYSISENUM - Enum of Balancethickness2Analysis
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=Balancethickness2AnalysisEnum()
+
+macro=StringToEnum('Balancethickness2Analysis');
Index: /issm/trunk-jpl/src/m/enum/Balancethickness2SolutionEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/Balancethickness2SolutionEnum.m	(revision 17831)
+++ /issm/trunk-jpl/src/m/enum/Balancethickness2SolutionEnum.m	(revision 17831)
@@ -0,0 +1,11 @@
+function macro=Balancethickness2SolutionEnum()
+%BALANCETHICKNESS2SOLUTIONENUM - Enum of Balancethickness2Solution
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=Balancethickness2SolutionEnum()
+
+macro=StringToEnum('Balancethickness2Solution');
Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py
===================================================================
--- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 17830)
+++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py	(revision 17831)
@@ -212,4 +212,5 @@
 def MeshVertexonbaseEnum(): return StringToEnum("MeshVertexonbase")[0]
 def MeshVertexonsurfaceEnum(): return StringToEnum("MeshVertexonsurface")[0]
+def MeshVertexonboundaryEnum(): return StringToEnum("MeshVertexonboundary")[0]
 def MeshXEnum(): return StringToEnum("MeshX")[0]
 def MeshYEnum(): return StringToEnum("MeshY")[0]
@@ -317,4 +318,6 @@
 def BalancethicknessAnalysisEnum(): return StringToEnum("BalancethicknessAnalysis")[0]
 def BalancethicknessSolutionEnum(): return StringToEnum("BalancethicknessSolution")[0]
+def Balancethickness2AnalysisEnum(): return StringToEnum("Balancethickness2Analysis")[0]
+def Balancethickness2SolutionEnum(): return StringToEnum("Balancethickness2Solution")[0]
 def BalancethicknessSoftAnalysisEnum(): return StringToEnum("BalancethicknessSoftAnalysis")[0]
 def BalancethicknessSoftSolutionEnum(): return StringToEnum("BalancethicknessSoftSolution")[0]
@@ -520,4 +523,5 @@
 def HydrologyWaterVxEnum(): return StringToEnum("HydrologyWaterVx")[0]
 def HydrologyWaterVyEnum(): return StringToEnum("HydrologyWaterVy")[0]
+def SigmaNNEnum(): return StringToEnum("SigmaNN")[0]
 def StressTensorEnum(): return StringToEnum("StressTensor")[0]
 def StressTensorxxEnum(): return StringToEnum("StressTensorxx")[0]
Index: /issm/trunk-jpl/src/m/enum/MeshVertexonboundaryEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/MeshVertexonboundaryEnum.m	(revision 17831)
+++ /issm/trunk-jpl/src/m/enum/MeshVertexonboundaryEnum.m	(revision 17831)
@@ -0,0 +1,11 @@
+function macro=MeshVertexonboundaryEnum()
+%MESHVERTEXONBOUNDARYENUM - Enum of MeshVertexonboundary
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=MeshVertexonboundaryEnum()
+
+macro=StringToEnum('MeshVertexonboundary');
Index: /issm/trunk-jpl/src/m/enum/SigmaNNEnum.m
===================================================================
--- /issm/trunk-jpl/src/m/enum/SigmaNNEnum.m	(revision 17831)
+++ /issm/trunk-jpl/src/m/enum/SigmaNNEnum.m	(revision 17831)
@@ -0,0 +1,11 @@
+function macro=SigmaNNEnum()
+%SIGMANNENUM - Enum of SigmaNN
+%
+%   WARNING: DO NOT MODIFY THIS FILE
+%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
+%            Please read src/c/shared/Enum/README for more information
+%
+%   Usage:
+%      macro=SigmaNNEnum()
+
+macro=StringToEnum('SigmaNN');