Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 25444)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 25445)
@@ -520,7 +520,5 @@
 		/*Start looping on Gaussian points*/
 		Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
-		for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-			gauss->GaussPoint(ig);
+		while(gauss->next()){
 			thickness_input->GetInputValue(&thickness,gauss);
 			calvingratex_input->GetInputValue(&calvingratex,gauss);
@@ -648,7 +646,5 @@
 		/*Start looping on Gaussian points*/
 		Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
-		for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-			gauss->GaussPoint(ig);
+		while(gauss->next()){
 			thickness_input->GetInputValue(&thickness,gauss);
 			calvingratex_input->GetInputValue(&calvingratex,gauss);
@@ -728,7 +724,5 @@
 	/* Start  looping on the number of gaussian points: */
 	gauss=new GaussPenta(0,1,2,2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		/*Compute strain rate viscosity and pressure: */
@@ -2050,7 +2044,5 @@
 	/*Start looping on Gaussian points*/
 	Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
 		vx_input->GetInputValue(&vx,gauss);
@@ -2172,7 +2164,5 @@
 	/*Start looping on Gaussian points*/
 	Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
 		vx_input->GetInputValue(&vx,gauss);
@@ -2289,6 +2279,6 @@
 			}
 			/*Integrate over edge*/
-			for(int ig=gauss[iv]->begin();ig<gauss[iv]->end();ig++){
-				gauss[iv]->GaussPoint(ig);
+			gauss[iv]->Reset();
+			while(gauss[iv]->next()){
 				penta->JacobianDeterminantLine(&Jdet,&xyz_list_line[0][0],gauss[iv]);
 				original_input->GetInputValue(&value,gauss[iv]);
@@ -3619,6 +3609,5 @@
 		for(int ig=0;ig<3;ig++){
 			GaussPenta* gauss=new GaussPenta(ig,3+ig,11);
-			for (int iv=gauss->begin();iv<gauss->end();iv++){
-				gauss->GaussPoint(iv);
+			while(gauss->next()){
 
 				/* Get the value we need*/
@@ -3839,7 +3828,5 @@
 	/*Start looping on Gaussian points*/
 	Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
 		calvingratex_input->GetInputValue(&calvingratex,gauss);
@@ -3966,7 +3953,5 @@
 	/*Start looping on Gaussian points*/
 	Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
 		calvingratex_input->GetInputValue(&calvingratex,gauss);
@@ -4015,7 +4000,5 @@
 	/* Start  looping on the number of gaussian points: */
 	gauss = this->NewGauss(point1,fraction1,fraction2,1-mainlyfloating,3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		this->JacobianDeterminantBase(&Jdet,&xyz_list[0][0],gauss);
 		floatingmelt_input->GetInputValue(&floatingmelt,gauss);
@@ -4060,7 +4043,5 @@
 	/* Start  looping on the number of gaussian points: */
 	gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		this->JacobianDeterminantBase(&Jdet,&xyz_list[0][0],gauss);
 		groundedmelt_input->GetInputValue(&groundedmelt,gauss);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 25444)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 25445)
@@ -672,7 +672,5 @@
 		/*Start looping on Gaussian points*/
 		Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
-		for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-			gauss->GaussPoint(ig);
+		while(gauss->next()){
 			thickness_input->GetInputValue(&thickness,gauss);
 			calvingratex_input->GetInputValue(&calvingratex,gauss);
@@ -815,7 +813,5 @@
 		/*Start looping on Gaussian points*/
 		Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
-		for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-			gauss->GaussPoint(ig);
+		while(gauss->next()){
 			thickness_input->GetInputValue(&thickness,gauss);
 			calvingratex_input->GetInputValue(&calvingratex,gauss);
@@ -2645,7 +2641,5 @@
 	/*Start looping on Gaussian points*/
 	Gauss* gauss=this->NewGauss(xyz_list,xyz_front,3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
 		vx_input->GetInputValue(&vx,gauss);
@@ -2774,7 +2768,5 @@
 	/*Start looping on Gaussian points*/
 	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
 		vx_input->GetInputValue(&vx,gauss);
@@ -2901,7 +2893,5 @@
 	/*Start looping on Gaussian points*/
 	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
 		vx_input->GetInputValue(&vx,gauss);
@@ -3367,6 +3357,5 @@
 	volume=0;
 
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		this->JacobianDeterminant(&Jdet,xyz_list,gauss);
@@ -3466,7 +3455,5 @@
 	/* Start  looping on the number of gaussian points: */
 	gauss=new GaussTria(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		/* Get Jacobian determinant: */
@@ -3505,7 +3492,5 @@
 	/* Start  looping on the number of gaussian points: */
 	gauss=new GaussTria(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 
 		/* Get Jacobian determinant: */
@@ -4381,7 +4366,5 @@
 	/*Start looping on Gaussian points*/
 	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
 		calvingratex_input->GetInputValue(&calvingratex,gauss);
@@ -4516,7 +4499,5 @@
 	/*Start looping on Gaussian points*/
 	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		thickness_input->GetInputValue(&thickness,gauss);
 		calvingratex_input->GetInputValue(&calvingratex,gauss);
@@ -4562,7 +4543,5 @@
 	/* Start  looping on the number of gaussian points: */
 	gauss = this->NewGauss(point1,fraction1,fraction2,1-mainlyfloating,3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		this->JacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
 		floatingmelt_input->GetInputValue(&floatingmelt,gauss);
@@ -4607,7 +4586,5 @@
 	/* Start  looping on the number of gaussian points: */
 	gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
+	while(gauss->next()){
 		this->JacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
 		groundedmelt_input->GetInputValue(&groundedmelt,gauss);
@@ -5770,7 +5747,6 @@
 		total_weight=0;
 		I=0;
-		for(int ig=gauss->begin();ig<gauss->end();ig++){
+		while(gauss->next()){
 			IssmDouble Ig=0;
-			gauss->GaussPoint(ig);
 			deltathickness_input->GetInputValue(&Ig,gauss);
 			I+=Ig*gauss->weight;