Index: issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/AdjointBalancethickness2Analysis.cpp	(revision 25439)
@@ -74,6 +74,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -175,6 +174,5 @@
 
 	Gauss* gauss=element->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -224,6 +222,5 @@
 
 	Gauss* gauss=element->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -278,6 +275,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -328,6 +324,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp	(revision 25439)
@@ -94,6 +94,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -288,6 +287,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		adjoint_input->GetInputDerivativeValue(&Dlambda[0],xyz_list,gauss);
@@ -335,6 +333,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		adjoint_input->GetInputDerivativeValue(&Dlambda[0],xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp	(revision 25439)
@@ -102,6 +102,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(5);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -179,6 +178,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(5);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -285,6 +283,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -397,6 +394,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGaussTop(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
@@ -636,6 +632,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGaussTop(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
@@ -889,6 +884,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -1228,6 +1222,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -1321,6 +1314,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		thickness_input->GetInputValue(&thickness,gauss);
@@ -1413,6 +1405,5 @@
 	/* 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);
+	while(gauss->next()){
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
@@ -1481,6 +1472,5 @@
 	/* 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);
+	while(gauss->next()){
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctionsP1(basis,gauss);
@@ -1542,6 +1532,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		thickness_input->GetInputValue(&thickness,gauss);
@@ -1634,6 +1623,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		thickness_input->GetInputValue(&thickness,gauss);
@@ -1729,6 +1717,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -1807,6 +1794,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGaussBase(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		adjointx_input->GetInputValue(&lambda, gauss);
@@ -1903,6 +1889,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGaussBase(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		adjointx_input->GetInputValue(&lambda, gauss);
@@ -2008,6 +1993,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		adjointx_input->GetInputValue(&lambda, gauss);
@@ -2084,6 +2068,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGaussBase(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		adjointx_input->GetInputValue(&lambda, gauss);
@@ -2178,6 +2161,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGaussBase(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		adjointx_input->GetInputValue(&lambda, gauss);
@@ -2283,6 +2265,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		adjointx_input->GetInputValue(&lambda, gauss);
@@ -2387,6 +2368,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		thickness_input->GetInputValue(&thickness,gauss);
Index: issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp	(revision 25439)
@@ -90,6 +90,5 @@
 	/* 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);
+	while(gauss->next()){
 		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -148,6 +147,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp	(revision 25439)
@@ -186,6 +186,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -284,6 +283,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -361,6 +359,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -402,6 +399,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -468,6 +464,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		thickness_input->GetInputValue(&thickness, gauss);
Index: issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp	(revision 25439)
@@ -106,6 +106,5 @@
 	/*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);
+	while(gauss->next()){
 
 		H_input->GetInputValue(&thickness,gauss);
@@ -200,6 +199,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp	(revision 25439)
@@ -472,6 +472,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -655,6 +654,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -749,6 +747,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -795,7 +792,5 @@
 	/* 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);
-
+	while(gauss->next()){
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctions(basis,gauss);
Index: issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp	(revision 25439)
@@ -70,6 +70,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -118,6 +117,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 25439)
@@ -138,6 +138,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=workelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		workelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp	(revision 25439)
@@ -69,6 +69,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp	(revision 25439)
@@ -69,6 +69,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp	(revision 25439)
@@ -183,6 +183,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -320,6 +319,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp	(revision 25439)
@@ -159,6 +159,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=topelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -298,6 +297,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=topelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp	(revision 25439)
@@ -103,6 +103,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp	(revision 25439)
@@ -215,6 +215,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss = basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss       ->GaussPoint(ig);
+	while(gauss->next()){
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
@@ -339,6 +338,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss = basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		basalelement ->NodalFunctions(basis,gauss);
Index: issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp	(revision 25439)
@@ -245,7 +245,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-
-	for(int ig=gauss -> begin();ig<gauss->end();ig++){
-		gauss       ->GaussPoint(ig);
+	while(gauss->next()){
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
@@ -384,7 +382,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		basalelement->NodalFunctions(basis,gauss);
Index: issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp	(revision 25439)
@@ -177,6 +177,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -255,6 +254,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp	(revision 25439)
@@ -102,6 +102,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -172,6 +171,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/L2ProjectionEPLAnalysis.cpp	(revision 25439)
@@ -132,6 +132,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -203,6 +202,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 25439)
@@ -308,6 +308,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -551,5 +550,5 @@
 
 	/*Intermediaries */
-	int i, ig, domaintype;
+	int domaintype;
 	IssmDouble  Jdet,dt;
 	IssmDouble  lsf;
@@ -573,7 +572,5 @@
 		/* Start  looping on the number of gaussian points: */
 		Gauss* gauss=basalelement->NewGauss(2);
-		for(ig=gauss->begin();ig<gauss->end();ig++){
-			gauss->GaussPoint(ig);
-
+		while(gauss->next()){
 			basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
 			basalelement->NodalFunctions(basis,gauss);
@@ -581,5 +578,5 @@
 			/* old function value */
 			levelset_input->GetInputValue(&lsf,gauss);
-			for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*lsf*basis[i];
+			for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*lsf*basis[i];
 		}
 
Index: issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 25439)
@@ -328,6 +328,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -517,6 +516,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -642,6 +640,5 @@
 
 	/* Start  looping on the number of gaussian points: */
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -748,6 +745,5 @@
 
 	/* Start  looping on the number of gaussian points: */
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -953,6 +949,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -1006,6 +1001,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -1079,6 +1073,5 @@
 
 	/* Start  looping on the number of gaussian points: */
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/MeltingAnalysis.cpp	(revision 25439)
@@ -102,6 +102,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=basalelement->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/SmoothAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/SmoothAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/SmoothAnalysis.cpp	(revision 25439)
@@ -74,6 +74,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -155,6 +154,5 @@
 	/* 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);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp	(revision 25439)
@@ -444,6 +444,5 @@
 
 		Gauss* gauss=element->NewGaussLine(nodedown,nodeup,3);
-		for(int ig=gauss->begin();ig<gauss->end();ig++){
-			gauss->GaussPoint(ig);
+		while(gauss->next()){
 
 			B_input->GetInputValue(&B,gauss);
Index: issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 25439)
@@ -354,6 +354,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGaussBase(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
@@ -412,6 +411,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -564,6 +562,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss   = element->NewGaussBase(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
@@ -623,6 +620,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGaussBase(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
@@ -683,6 +679,5 @@
 	/* Start  looping on the number of gaussian points: */
 	Gauss* gauss=element->NewGauss(4);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
Index: issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp	(revision 25438)
+++ issm/trunk-jpl/src/c/analyses/UzawaPressureAnalysis.cpp	(revision 25439)
@@ -87,7 +87,5 @@
 
 	//Gauss* gauss = element->NewGauss(5);
-	//for(int ig=gauss->begin();ig<gauss->end();ig++){
-	//	gauss->GaussPoint(ig);
-
+	//while(gauss->next()){
 	//	element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 	//	this->GetM(M,element,gauss);
@@ -136,6 +134,5 @@
 
 	Gauss* gauss = element->NewGauss(5);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -251,6 +248,5 @@
 		delete gauss;
 		gauss = element->NewGaussBase(3);
-		for(int ig=gauss->begin();ig<gauss->end();ig++){
-			gauss->GaussPoint(ig);
+		while(gauss->next()){
 
 			/*Compute Jlambda*/
