Index: issm/trunk/src/c/classes/FemModel.h
===================================================================
--- issm/trunk/src/c/classes/FemModel.h	(revision 21729)
+++ issm/trunk/src/c/classes/FemModel.h	(revision 22758)
@@ -20,6 +20,9 @@
 class Profiler;
 class Elements;
-#ifdef _HAVE_NEOPZ_
+#if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_)
 #include "./AdaptiveMeshRefinement.h"
+#endif
+#if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_)
+#include "./AmrBamg.h"
 #endif
 /*}}}*/
@@ -47,6 +50,11 @@
 		Vertices    *vertices;             //one set of vertices
 
-		#ifdef _HAVE_NEOPZ_
+		//FIXME: do we want only one class and have virtual functions? or keep 2 classes, at least rename AdaptiveMeshRefinement -> AmrNeopz
+		#if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_)
 		AdaptiveMeshRefinement *amr;		  //adaptive mesh refinement object. It keeps coarse mesh and execute refinement process
+		#endif
+
+		#if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_)
+		AmrBamg *amrbamg; //adaptive mesh refinement object. It keeps coarse mesh and execute refinement process
 		#endif
 
@@ -61,4 +69,5 @@
 		FemModel* copy();
 		void Echo();
+		int  GetElementsWidth(){return 3;};//just tria elements in this first version
 		void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, char* restartfilename, const int solution_type,bool trace,IssmPDouble* X=NULL);
 		void InitFromFids(char* rootpath, FILE* IOMODEL, FILE* toolkitsoptionsfid, int in_solution_type, bool trace, IssmPDouble* X=NULL);
@@ -73,5 +82,5 @@
 		/*Modules*/ 
 		void BalancethicknessMisfitx(IssmDouble* pV);
-		void CalvingRateDevx();
+		void CalvingRateVonmisesx();
 		void CalvingRateLevermannx();
 		void DeviatoricStressx();
@@ -79,10 +88,10 @@
 		void ElementOperationx(void (Element::*function)(void));
 		void ElementResponsex(IssmDouble* presponse,int response_enum);
-		void FloatingAreax(IssmDouble* pV);
+		void FloatingAreax(IssmDouble* pV, bool scaled);
 		void GetInputLocalMinMaxOnNodesx(IssmDouble** pmin,IssmDouble** pmax,IssmDouble* ug);
-		void GroundedAreax(IssmDouble* pV);
-		void IceMassx(IssmDouble* pV);
-		void IceVolumex(IssmDouble* pV);
-		void IceVolumeAboveFloatationx(IssmDouble* pV);
+		void GroundedAreax(IssmDouble* pV, bool scaled);
+		void IceMassx(IssmDouble* pV, bool scaled);
+		void IceVolumex(IssmDouble* pV, bool scaled);
+		void IceVolumeAboveFloatationx(IssmDouble* pV, bool scaled);
 		void MassFluxx(IssmDouble* presponse);
 		void MaxAbsVxx(IssmDouble* presponse);
@@ -98,11 +107,13 @@
 		void MinVyx(IssmDouble* presponse);
 		void MinVzx(IssmDouble* presponse);
+		void DistanceToFieldValue(int fieldenum,IssmDouble fieldvalue,int distanceenum);
 		void ResetLevelset();
 		void StrainRateparallelx();
 		void StrainRateperpendicularx();
+		void StrainRateeffectivex();
 		void StressIntensityFactorx();
-		void TotalFloatingBmbx(IssmDouble* pFbmb);
-		void TotalGroundedBmbx(IssmDouble* pGbmb);
-		void TotalSmbx(IssmDouble* pSmb);
+		void TotalFloatingBmbx(IssmDouble* pFbmb, bool scaled);
+		void TotalGroundedBmbx(IssmDouble* pGbmb, bool scaled);
+		void TotalSmbx(IssmDouble* pSmb, bool scaled);
 		#ifdef  _HAVE_DAKOTA_
 		void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);
@@ -116,4 +127,6 @@
 		void Responsex(IssmDouble* presponse,const char* response_descriptor);
 		void SurfaceAbsMisfitx( IssmDouble* pJ);
+		void OmegaAbsGradientx( IssmDouble* pJ);
+		void EtaDiffx( IssmDouble* pJ);
 		void ThicknessAbsGradientx( IssmDouble* pJ);
 		void ThicknessPositivex(IssmDouble* pJ);
@@ -122,5 +135,5 @@
 		#endif
 		#ifdef _HAVE_ESA_
-		void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* xx, IssmDouble* yy); 
+		void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy); 
 		void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); 
 		#endif
@@ -128,5 +141,5 @@
 		void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
 		void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution);
-		void SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
+		void SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
 		void SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); 
 		IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg);
@@ -137,4 +150,7 @@
 		void UpdateConstraintsExtrudeFromTopx();
 		void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count);
+		void InitTransientOutputx(int* input_enum, int numoutputs);
+		void StackTransientOutputx(int* input_enum,int* stackedinput_enum, IssmDouble hydrotime, int numoutputs);
+		void AverageTransientOutputx(int* input_enum,int* averagedinput_enum,IssmDouble hydrotime,int numoutputs);
 		void UpdateConstraintsx(void);
 		int  UpdateVertexPositionsx(void);
@@ -146,13 +162,11 @@
 		#endif
 
-		#ifdef _HAVE_NEOPZ_
-		/*Adaptive mesh refinement methods*/
-		void InitializeAdaptiveRefinement(void);
+		/*AMR*/
+		#if !defined(_HAVE_ADOLC_)
 		void ReMesh(void);
 		void BedrockFromMismipPlus(void);
 		void AdjustBaseThicknessAndMask(void);
 		void GetMesh(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist);
-		int GetElementsWidth(){return 3;};//just tria elements in this first version
-		void ExecuteRefinement(int &numberofvertices,int &numberofelements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist);
+		void GetMeshOnPartition(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist,int** psidtoindex);
 		void GetGroundediceLevelSet(IssmDouble** pmasklevelset);
 		void CreateVertices(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,int* my_vertices,IssmDouble* newx,IssmDouble* newy,IssmDouble* newz,Vertices* vertices);
@@ -160,13 +174,35 @@
 		void CreateMaterials(int newnumberofelements,bool* my_elements,Materials* materials);
 		void CreateNodes(int newnumberofvertices,int* my_vertices,int nodecounter,int analysis_enum,Nodes* nodes);
-		void CreateConstraints(int newnumberofvertices,int newnumberofelements,int nodecounter,int constraintcounter,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints);
+		void CreateConstraints(Vertices* newfemmodel_vertices,int nodecounter,int constraintcounter,int analysis_enum,Constraints* newfemmodel_constraints);
+		void GetInputs(int* pnumP0inputs,IssmDouble** pP0inputs,int** pP0input_enums,int** pP0input_interp,int* pnumP1inputs,IssmDouble** pP1inputs,int** pP1input_enums,int** pP1input_interp);
 		void InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements);
 		void UpdateElements(int newnumberofelements,int* newelementslist,bool* my_elements,int nodecounter,int analysis_counter,Elements* newelements);
 		void ElementsAndVerticesPartitioning(int& newnumberofvertices,int& newnumberofelements,int& elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices);
 		void WriteMeshInResults(void);
+		void WriteErrorEstimatorsInResults(void);
+		void SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy); //nodal values, just for SSA-P1: TauXX, TauYY, TauXY
+		void ZZErrorEstimator(IssmDouble** pelementerror);
+		void SmoothedGradThickness(IssmDouble** pdHdx,IssmDouble** pdHdy);
+		void ThicknessZZErrorEstimator(IssmDouble** pelementerror);
+		void MeanGroundedIceLevelSet(IssmDouble** pmasklevelset);
+		void GetElementCenterCoordinates(IssmDouble** pxc,IssmDouble** pyc);
+		void GetZeroLevelSetPoints(IssmDouble** pzerolevelset_points,int &numberofpoints,int levelset_type);
+		#endif
+
+		#if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_)
+		void ReMeshBamg(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist);
+		void InitializeAdaptiveRefinementBamg(void);
+		void GethmaxVerticesFromZeroLevelSetDistance(IssmDouble* hmaxvertices,int levelset_type);
+		void GethmaxVerticesFromEstimators(IssmDouble* hmaxvertices,int errorestimator_type);
+		void GetVerticeDistanceToZeroLevelSet(IssmDouble** pverticedistance,int leveset_type);
+		#endif
+
+		#if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_)
+		void ReMeshNeopz(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist);
+		void InitializeAdaptiveRefinementNeopz(void);
+		void GetElementDistanceToZeroLevelSet(IssmDouble** pelementdistance,int levelset_type);
 		void SetRefPatterns(void);
 		#endif
 };
-		
 
 #endif
